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ABSTRACT 


A numerical technique with simultaneous time integration 
of a meshed grid system is proposed, in which the fine mesh 
region is able to move within the coarse grid. The interface 
boundary conditions employed are shown analytically to be the 
only stable specification of those tested for a simple linear 
case. Numerical solutions using the meshed system are con- 
pared with those from uniform coarse and fine grids. One im- 
portant criterion is that the solution within the fine mesh 
region of the meshed grid must have nearly the same accuracy 
as in a system which uses a fine mesh everywhere. The tech- 
nique is applied to a two-dimensional (y,p) ten-level primi- 
tive equation model. Behavior of the meshed model is exam- 
ined in experiments in which a small scale heat source is 
imbedded within an undisturbed zonal flow pattern. Energy 
boundary fluxes across the interface in the meshed system 
are shown to compare favorably with those in a uniform fine 
mesh grid. Three-dimensional (x,y,p) versions of the primi- 
tive equation model with single and multiple grid reductions 
are also examined. Simulated small scale disturbances are 
used to compare the meshed model solutions with those obtained 
on a uniform grid. Future applications of the meshing tech- 
nique are suggested, such as in tropical storm research and 


general circulation models of the atmosphere and ocean. 
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I. INTRODUCTION 


Use of nested finite grids in numerical applications has 
been studied by many authors in various disciplines. Although 
the approach may vary in the treatment of this subject, the 
ultimate goal is the same, namely higher resolution in speci- 
fied areas of the grid network. The geophysical applications 
,of such grid systems are many and varied. In atmospheric 
general circulation models, for example, higher resolution 
may be desirable over continental areas where orographic 
effects are important. Ocean models, on the other hand, 
often have continental boundary regions which require more 
careful consideration compared to the large scale flow. At- 
mospheric prediction models which are generally unable to 
adequately resolve such phenomena as fronts and tropical 
storms could incorporate fine mesh grids in these regions. 
There are several basic methods of approaching the limited 
area forecasting problem. The first method is to design a 
model in which a large enough geographic area is covered to 
permit boundary condition errors to be neglected. Except 
for the initialization, there can be no effect on the in- 
terior forecast by changes in the planetary scale flow with 
this technique. This, of course, limits the usefulness of 
such specifications to very short range forecasts. In addi- 
tion, this method wastes both computer time and memory 
space. A second technique is to design a nested grid system 


in which the grid scale abruptly changes, and specifying the 
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required boundary values for the limited area grid from the 
larger scale forecast, or in some way assigning them from 
the interior at outflow points. It will be shown later that 
if equations on the grids are simultaneously integrated, 
there is a particular method of specifying boundary values 
at the interface that is computationally stable and produces 
a superior forecast. 

There are, however, two other approaches to the problem 
of increased grid resolution in certain areas. The first of 
these involves the mathematical transform of a non-orthogonal, 
unequal spatial increment grid into a constant interval, 
orthogonal grid. Anthes (1970) has reported such a technique 
in which the non-orthogonal grid to be transformed had the 
property that the mesh size in both north-south (y), and east- 
west (x) directions varied smoothly from a minimum value at 
the center of the domain to a maximum value at the boundary. 
The equations were then transformed in terms of new indepen- 
dent variables on an orthogonal mesh, solved in the new 
system, and the results then transformed back to the origi- 
nal domain. Anthes showed that such a variable grid may be 
used in certain physical problems to gain resolution at the 
center of the grid over a limited time range. The major dis- 
advantage of this type of formulation is the inability of 
the region of maximum resolution to move without moving the 
entire grid system. The small scale system of interest 
could well move a substantial distance from the high resolu- 


tion area over the forecast period. Two computational 





difficulties encountered with an analytic transformation 
are the loss of second order accuracy for re space 
differencing, and the addition of non-linear terms result- 
ing from the non-orthogonality of the transformed grid 
points in real space. Finally, one can retain orthogonality 
of a smoothly varying grid in one space dimension such as 
in axisymmetric hurricane models (Estoque 1962, Kuo 1965, 
Yamasaki 1968). Sundqvist and Veronis (1970) have shown 
that second order accuracy can be retained in a variable 
grid if the rate of variation is properly chosen. Since 
this technique does not apply to two-dimensional horizontal 
grids, it can only find application in a rather limited 
number of geophysical applications. 

The systematic approach to the design of a nested grid 
System should begin with simple one-dimensional experiments 
with results which are applicable to more general systems, 
as opposed to the smoothly varying grid which is only 
applicable in one space dimension. Hence, one must evaluate 
the features which are desired in the nested system in 
general, and apply them in a simple case. The main problem 
is obvious, namely, the interface boundary condition of the 
nested system. In the first section of this report, it shall 
be shown analytically and with numerical integrations of a 
Simple system, that only one of the types of boundary condi- 
tions tested can be specified when centered space and time 
differencing are used. This specification is computationally 


stable and leads to an accurate numerical solution. One- 





dimensional numerical experiments show that many specifica- 
tions which are computationally stable in uniform grid net- 
works are extremely unstable in the meshed system. In the 
design of the system, one must also insure that the fine 
mesh grid (FMG) is able to move within the coarse mesh grid 
(CMG). This is necessary because it is desirable to have 
the smallest FMG possible to conserve computer time and 
memory space. The ultimate objective is to produce a fore- 
cast as accurate as that obtained on a grid which is fine 
mesh everywhere, but at a fraction of the computer cost. 
Thus, several one-dimensional tests, both linear and nonlin- 
ear, were performed to demonstrate this capability of the 
meshed system. Phase and amplitude errors are compared to 
those in uniform grid systems, and the meshed system is 
shown to compare favorably to the uniform fine grid. 

A two-dimensional (y,p) primitive equation (PE) model 
was used in subsequent two-dimensional tests of the meshing 
technique. One of the primary objectives here was to study 
the boundary fluxes of kinetic and thermal energy through 
the interface. These fluxes could then be compared to those 
across the same region in a uniform fine grid. In these 
experiments, a localized heat source simulating a band of 
convective elements in an initially undisturbed flow was 
employed. Comparison of forecasts of various parameters were 
made between uniform coarse grids, uniform fine grids and the 
meshed grid model. In each case, the meshed system produced 


forecasts comparable to the uniform fine grid in the FMG 





region, and superior to the uniform coarse grid outside this 
region. 

Final tests of the meshing technique were made utilizing 
a four-level primitive دھگ۔۔‎ model, suitable for conver- 
sion to possible operational uses with minimal changes. This 
model had two basic versions, one in which the grid scale 
changed abruptly from the coarse grid to the fine, and 
another in which the scale still changed abruptly, but in two 
steps. Thus, in the latter version, there were actually 
three different grid scales. Comparisons of uniform fine 
grid solutions were not possible for any of the three-dimen- 
sional cases because of computer memory and time limitations. 
Analytic initial states which simulated tropical storm and 
frontal scale disturbances were used in the comparisons of 
various aspects of the forecast fields. As in the two-dimen- 
sional model, realistic fluxes of kinetic and thermal energy 
were computed at the interface. Movement of the FMG while 
following the simulated disturbance had no apparent deleter- 
ious effects on the performance of the meshed system. Sugges- 
ted future uses of the three-dimensional model are discussed 


as well as problems involving initialization of the FMG. 
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II. REVIEW OF MESHING TECHNIQUES 


One of the earliest applications of the concept of abrupt- 
ly changing grid scales in a prediction model was made by 
Birchfield (1960). Utilizing the simple barotropic vorticity 
model, he showed improved trajectory forecasts of several 
tropical storms could be attained with a meshed grid system. 
The grid scale of the fine mesh grid in Birchfield's model 
was one-half that of the coarse mesh grid, and the fine mesh 
Was initialized by interpolation from the CMG. One of the 
unique features of Birchfield's model was the simultaneous 
time integration of both grid scales. In the marching pro- 
cess, the vorticity equation was solved by relaxation tech- 
niques over the entire CMG. During this step, the presence 
of the FMG was ignored. After the CMG time step was complete, 
spatial and temporal interpolation of the tendencies at co- 
incident points yielded boundary values for the subsequent 
FMG time steps. In each grid, the time increment was such 
that the stability criterion of Courant et al (1928) was 
satisfied. Thus, four time steps were taken in the FMG for 
everyone in the CMG. Birchfield used one-fourth the time 
increment in the FMG because of the higher wind speeds in 
the storms. Interaction between scales was attained in this 
model by re-initializing the CMG after each time step, hold- 
ing the FMG predicted value fixed at interior coincident 
points. This model, although one of the earliest meshed 


models, incorporated two very important features, namely 
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simultaneous time integration and a movable FMG. It will be 
shown later that simultaneous time integration on both grids 
is necessary for a stable, accurate forecast. Phillips 
(1972) also argued in favor of the "strategy" of simultaneous 
time integration. 

Little concerning grid nesting techniques appeared in the 
literature until Howcroft's (1966) description of the proposed 
Air Weather Service's (AWS) fine mesh, barotropic vorticity 
model. The FMG spatial increment in the AWS model was again 
one-half that of the CMG. Initial data for the FMG were ob- 
tained by interpolation of geopotential heights at the time 
of Howcroft's report, however he noted that a fine grid anal- 
ysis program was under development. The integration made use 
of a 24-hour forecast on a hemispheric coarse grid, retaining 
the appropriate coincident point values for use as boundary 
conditions in the subsequent FMG forecast. The fine mesh 
forecast was then carried out, updating the boundary values 
from the previous hemispheric forecast. This became a rather 
popular way of specifying boundary values for nested models, 
but it will be shown later that the method is generally un- 
Satisfactory, even in linear models. 

The need for detailed forecasts in active regions is re- 
viewed in some detail by Bushby and Timpson (1967). They 
noted the important point that if fine mesh integrations are 
started from smooth realistic patterns of motion, temperature 
and humidity, obtained from large scale analysis, one may 


í 
expect smaller scale detail to develop if the problem is 
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properly formulated. This is crucial for fine mesh models, 
for if it is not true, models with grid increments smaller 
than the disfance between sounding stations have little mean- 
ing for operational forecasting. Bushby and Timpson's model 
was fine mesh everywhere, with the lateral boundary conditions 
specified from the interior of the forecast region. Such 
specifications will later be shown computationally unstable 

if not very carefully treated, and even in the stable cases 
may develop severe resonance. In a later report, Bushby(1968) 
noted that specification from the interior was in fact un- 
satisfactory, and that temporally constant lateral boundaries 
were adopted with an enhanced diffusion coefficient near the 
boundaries. 

The use of independent large scale (hemispheric) fore- 
casts to specify FMG boundaries was also used by Hill (1963) 
in the so-called "telescope" model. A numerical integration 
over an octagonal hemispheric grid was first carried out. 
Another grid with the same shape and number of points, but 
with half the grid spacing, was placed within the original 
grid, so that about one-quarter of the grid points coincided. 
The forecasts were generated separately, except for the 
boundary values for the FMG, which were obtained by spatial 
and temporal interpolation from the e eric forecast. 

The process may be repeated by placing a third grid within 
the second grid and so on. A later report concerning Hill's 
method by Shapiro and O'Brien (1970) proposed that a type of 


upstream difference be taken at outflow points, with the . 
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previously discussed specification at the inflow points re- 
tained. This technique, originally suggested by Charney 
(1960), will be shown to be satisfactory, under certain con- 
ditions. 

An important natural advantage of a fine mesh grid is 
the opportunity to portray the effects of the boundary layer 
in much greater detail. Horizontal variations due to features 
such as orography and boundary fluxes of heat and moisture 
can be simulated in a more realistic fashion as the grid 
scale begins to approach the scale of the phenomena themselves. 
A common feature in many operational prediction and general 
circulation models is the crude treatment of orographic fea- 
tures, particularly extensive mountain ranges such as the 
Alps, which have a profound effect on weather systems. Ros- 
seau (1968), has reported on a meshed grid system developed 
for detailed local forecasts in the region of the Alps and 
Pyrenees. At some predetermined time (T hours before veri- 
fication time) during the integration, the data of a hemis- 
pheric forecast in the region of interest was extracted from 
the forecast fields and interpolated to the FMG. This became 
the initial state for the regional forecast model over west- 
ern Europe, which had a grid increment one-tenth that of the 
hemispheric model. Very detailed orography was modeled to 
produce regional forecasts for France during the FMG inte- 
gration over the T remaining hours. Rosseau noted that ex- 
tracting boundary values from the subsequent independent 


hemispheric forecast led to severe computational noise, 
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which masked meteorological development. As a consequence, 
Rosseau formulated outflow boundary conditions based on the 
upstream differencing of the potential vorticity, as sugges- 
ted by Charney (1960). The major disadvantage of this spe- 
cification is the lack of interaction between the grid scales. 
Due to the greatly increased computer time and memory space 
required for such a fine grid, this integration was only for 
three hours for one data set, making evaluation difficult. 

Another approach to the limited area forecasting problem 
was that by Gerrity and McPherson (1969). In this barotropic 
primitive equation (PE) model, the integration was carried 
out on a uniform fine grid mainly covering the North Ameri- 
can continent. The grid spacing for this model was about 
190 km at 60N, or one-half that of the National Meteorologi- 
cal Center (NMC) hemispheric model. The boundary values were 
temporally constant, and Euler “backward" differencing was 
used throughout. Since the "semi-momentum" (Shuman and 
Hovermale 1968) space. differencing was used, as well as 
external smoothers, it is difficult to determine the effects 
of the temporally constant boundary values in this model. 
It is clear, however, that such a model is useful only for 
very short range forecasts in regions well away from the 
boundary. 

Wang and Halpern (1970) reported experiments with a 
barotropic PE model similar to that just described. The 
spatial differencing was the same, but the boundary condi- 


tions were obtained from an independent CMG forecast, and 


15 





leapfrog-trapezoidal (Kurihara 1965) time differencing was 
employed. External smoothers were also employed to remove 
the noise inherent in the boundary specification. The major 
disadvantage was absence of interaction between grid scales. 
A more sophisticated use of the barotropic primitive equa- 
tion model was made by Ookochi (1972), in which he employed 
simultaneous time integration on both scales. This model 
was solved on a staggered grid using Euler backward time 
differencing, and a spatial averaging technique to describe 
the boundary conditions at the interface between grid scales. 
Ookochi experimentally showed the method to be computation- 
ally stable over the period of integration, but the model 
does not include the type of gravity waves discussed by 
Phillips (1972). Another disadvantage was the lack of pro- 
vision for moving the fine grid region, resulting in distor- 
tion of phase speed as relatively short waves passed from 


one grid to the other. 
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III. THE MESHING TECHNIQUE AND BOUNDARY 


EONDERFONTEXPERIMENTS 


A. GENERAL CONSIDERATION 

The most obvious problem in a nested grid system is the 
treatment of the interface between the two scales. Mathema- 
tical treatment of a similar problem was given by Osher 
(1970) for solution of the heat diffusion equation on a grid 
with an abrupt spatial increment change. Appropriate differ- 
ence equations were shown to be maximum norm stable over the 
whole spatial domain, and particularly at the interface. 
Ciment (1971) examined the solution of the linear equation 


du _ 
> U 


(1) 


a 
le 


across a discontinuity in grid spatial increment and showed 
certain dissipative differencing schemes such as the Lax- 
Wendroff, are stable approximations to (1) over the whole 
domain. It will be shown that under certain conditions, the 
centered space and time differencing commonly used in met- 
eorological models is also a stable approximation. Treat- 
ment of the outflow boundary for (1) has been studied by 

many authors. Outflow boundaries are particularly trouble- 
some when centered spatial differencing is used due to pos- 
sible separation of the solution at alternate grid points. 
Platzman (1954) showed that, in general, the outflow boundary 
value should be externally specified, although there are cer- 
tain cases which allow a specification from the interior. 


Many possible outflow boundary specifications have been 
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examined in numerical tests by Nitta (1962), who, however, 
did not include a test in which the exact solution was spe- 
cified on the outflow boundary. It will be shown later that 
not only is this specification computationally stable, but 
yields an accurate numerical solution for low wave numbers. 
Since this specification is possible for only very simple 
cases, one may argue heuristically that specifying a numer- 
meat solution from a larger grid at the outflow point might 
also give good results. In the nested system, the problem 
is compounded, as outflow boundaries for both grid scales 


Must be considered. 


B. COMPUTATIONAL STABILITY CONSIDERATIONS 

Three basic grids were employed for the linear tests to 
demonstrate the meshing technique and compare various bound- 
ary specifications. The first grid was a uniform mesh con- 
trol which merely specified cyclic continuity of the periodic 
initial state. Numerical solutions of (1) on such a grid 
are thus free of boundary conditions, and if the stability 
criterion of Courant et al (1928) is met, the solution can- 
not become contaminated by computational noise, This grid 
was identical to the CMG grid in the meshed model, and 
direct comparisons of numerical solutions were made. Thus, 
the only factor which could possibly have caused differences 
was the presence of the FMG. The second grid system utilized 
was similar to that of Nitta (1962), wherein the exact solu- 
tion is specified at the inflow, while various outflow con- 


ditions may be tested in the numerical solution of (1). 
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This grid, which was thus quite different from the nested 
model which was cyclic, was very convenient in the outflow 
condition tests. The third grid was, of course, the meshed 


grid system, an example of which is depicted in Figure 1. 


j=1 2 3 45678910 N 12 13 14 





Fig. 1. Meshed grid for p = 2. The CMG is cyclic, and 
for this example, the FMG covers three CMG 
intervals. 

Note that the grid abruptly changes spatial increment at 

the points j = 4 and j = 10. An integer p is defined, 

which is the ratio of the spatial increment (Ax,) in the 

CMG to the increment (Ax, ) in the FMG. | 

An approach similar to that of Platzman (1954)was util- 
ized for the meshed model, to demonstrate computational 
stability of a certain interface boundary specification 
which will henceforth be called "mutually determined." The 


mutually determined boundary condition is specified in the 


following manner. In Figure 1, for example, the solution 


domain of the FMG is defined to be from j = 4 to j = 10, 
and the CMG domain from j = 1 to j = 6, and from j = 8 to 
j = 14. Thus, there is an overlap of the solution domains 


of one CMG grid increment. In the simultaneous time 
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integration, the solution at the points j = 4 and j = 10 in 
Figure 1 is determined by the CMG, while at j = 6 and j = 8, 
the solution is on the FMG. In this fashion, a dual initial- 
boundary problem may be solved with simultaneous time inte- 
gration. Because the solution from each grid is used to 
determine boundary conditions for the other, mutual inter- 
action can occur across the interface. There are, of course, 
other specifications which may be utilized at both ends of 
the FMG, but for the following tests, only the FMG outflow 
boundary was modified. On the inflow point, the mutually 
determined condition was retained. The main point of these 
tests was to show that the only computationally stable spe- 
cification of those tested for the outflow boundary of the 
FMG was the mutually determined condition. 

Approximating the spatial derivative of (1) with a cen- 
tered difference, but keeping time as a continuous variable, 
yields the following non-dimensional, "semi-discrete" dif- 


ference equation 





du 1 
یی ۳ئ تہ‎ Dun 
in the CMG, and 
ei, l 
ee) (209 
in the FMG, where t=- غال‎ . The external boundaries of the 
2 


CMG are cyclic, so that any contamination of a numerical 
solution must arise at the interface. The point j = 10 is 
thus the outflow boundary for the FMG if it is assumed that 


U > 0 in (1), and various boundary conditions were specified 
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Bere. At the outflow point for the FMG (j = 10 in Figure 1) 


the mutually determined condition is 





z ul] 7 "8? 


in the CMG, where Us is calculated on the FMG, and 


in the FMG. The mutually determined condition at the inflow 


point j = 4 is 
OTT 
oT 2 6 3 
in the CMG and 
dus _1 oT 
IT 4 6 4 


in the FMG. The system of difference equations over the 
whole spatial domain for this simple linear example may be 
expressed by the matrix equation 


9 
مو‎ AN 3) 


“ 
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where the first and last rows reflect ene ey Clic, continuity 
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in the CMG. The factor .> appears in the equations for the 
FMG because the space increment is one-half that of the CMG 
for this example. The fourth and tenth rows indicate the 
CMG difference equations at j = 4 and j = 10 respectively, 
while the adjacent interior rows represent the FMG differ- 
ence equations at j = 5 and j = 9. 

The matrix equation (3) is a matrix differential equa- 
tion of the first order, with the fundamental solution exp 
At. If the eigenvalues of A are ۸۰ ۰۰-۰٠٦۷٠٥۸... ) the 
fundamental modes of u are expA,t, eXpAnt,... These semi- 
discrete modes will then amplify or damp if the real parts 
of the A, are positive or negative, respectively. If the 
eigenvalues have no real parts, the semi-discrete modes are 
neutral. Platzman (1954) has shown that if any of the semi- 
crete modes either amplify or damp, then the modes of 
the fully discrete system 


n+l n-l n n 
u. = u. + AT (u. - u. (4a) 
j j ( F1 7 

in the CMG, and 


T 
untl s anl g کے‎ 


n n 
. 4b 
J J DD 


او ا کت 

in the FMG, will amplify exponentially. This was defined 

as the "secondary" stability criterion by Platzman (1954), 
as opposed to the "primary" criterion developed by Charney 
ه‎ 41 (1950), which requires that the magnitude of the maxi- 
mum eigenvalue of the fully discrete amplification matrix 

be equal to, or less than unity. The eigenvalues of A may 


be obtained numerically, and for this case, all the 
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eigenvalues of A were pure imaginary numbers. Thus, in the 
fully discrete integration of (1) with centered space and 
time differencing, one may expect no exponential growth due 
to the boundary specification. To further illustrate this 
Stability for the meshed system, the fully discrete amplifi- 
cation matrix may be obtained in the following manner. In 


(4), define a new variable "> such that 


The fully discrete system may now be written 


where for the example in Figure 1 


W col (W1 » War War...) Uy» Unser.) 


and 


us 
N 





where I is the identity matrix and 0 is the zero matrix. 
Again, the eigenvalues for B were numerically determined, 
with the largest magnitude equal to 1.0000. Thus, the 
stability criterion derived from the semi-discrete system 
is also demonstrated for the fully-discrete case. 
C. ANALYTIC AND NUMERICAL TESTS WITH OUTFLOW BOUNDARY 
CONDITIONS 
The semi-discrete stability analysis was then applied 
to several other specifications in the meshed model, as well 


as in a uniform grid system. For specifications involving 
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more than one time level, such as the upstream differencing, 
and for external specifications, numerical tests using cen- 
tered space and time differencing were performed. Table 1 
Summarizes the results of these tests. In the listing of 
specifications in Table 1, the outflow point index is j, 
while n = 0,1,2... is the time level for the integration of 
the fully discrete system (4). For the numerical tests, the 
simple initial state was a sinusoid of wave number 1 on the 
one-dimensional grid of unit non-dimensional length. In 
Table 1, t' is defined in the usual manner where x = Lx', 
and t = L/Ut', if the scale of u is U in (1). The one-di- 
mensional grid utilized a CMG spatial increment of 0.01 non- 
dimensional units and an FMG increment of 0.005 non-dimen- 
sional units. In the uniform grid, the spatial increment 
was the same as that in the CMG, and a time increment of 
0.0025 non-dimensional units was used for all cases, regard- 
less of the spatial increment of the grid. Later in this 
section, we examine the behavior of the "stable" specifica- 
tions with initial states involving higher wave numbers. 
The various semi-~discrete amplification matrices for the 
cases in Table 1 may be found in Appendix A. There is not 
a one-to-one correspondence between the magnitude of the 
largest real eigenvalue in Table 1 and the time at which 
exponential growth of the solution initiated because of the 
large number of computational modes for each solution. 
Specification 1 (mutually determined) is, of course, not 


applicable to the uniform grid, and as seen earlier, no 
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Table 1 


Summary of outflow Boundar tests 











1 Largest alu Critical a 
Outflow Specification Uniform Grád Meshed Grid | Uniform Grid, Meŝhed Crid 
Mutually determined 0.0000 stable 
Independent CMG stable 
sin 2T[(3Ax)-(nm+1)A7] stable 
0.0000 0.1269 stable 
1.0000 0.2423 fe 8, 
0.0000 OMS 32 stable 
1.5000 0.3565 2,8 
stable 
i سے‎ 
See text for explanation. 
2 


Largest real part of the maximum eigenvalue of the semi-discrete amplification 
matrix when such matrix is defined. 


time at which exponential growth began in the fully discrete‏ کر ٹر کت 


numerical integration. 
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instabilities arise from this interface condition. The nu- 

merical solution for this specification after 4000 time steps 
(t' = 10) is depicted in the upper portion of Figure 2, com- 
pared with the exact solution. The location of the station- 


ary FMG is denoted by the vertical dashed lines in Figure 2. 


AMPLITUDE 
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10 


Fig. 2. Numerical (solid) and exact (dashed) solutions at 
t' = 10 for specification l1l (upper) and specifica- 
tion 2 (lower). See text and Table 1 for details. 

The solution exhibited no computational noise, amplitude 

was conserved to within 0.01 percent, and only a small phase 


error occurred for this relatively long wave. This may be 


compared to the solution portrayed in the lower portion of 
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Figure 2, which results when an independently integrated CMG 
(specification 2 in Table 1) is used to specify the bound- 
ary values of FMG. This non-simultaneous solution corres- 
ponds to the approach of solving the equations on a coarse 
grid first, independent of the fine The appropriate 
coincident point values are retained for use as boundary 
values in the subsequent FMG integration. As the solutions 
on the different grids diverge, slight inconsistencies in 
the FMG boundary values quickly fill the forecast region 
with computational noise. Because this is essentially an 
external specification, no exponential growth of the compu- 
tational modes occurred, in fact, amplitude was conserved 

to within two percent. As seen in Figure 2, however, the 
solution is so contaminated that even the use of heavy 
smoothing would not render it usable. In a nonlinear model, 
this computational noise would lead to en. and pos- 
sibly to nonlinear computational instability. Phillips 
(1972) has shown, using the method of characteristics, that 
the independently integrated CMG specification can probably 
not be applied in a general multi-level model. In simple 
barotropic models, which do not allow internal gravity waves, 
phase speed magnitude of external gravity waves is always 
greater than the advective speed, but in multi-level primi- 
tive equation models, values of internal gravity wave pro- 
pagation speed will be encountered for short vertical wave- 
lengths which are less than the advective speed. For this 


case, Phillips showed that no external specification should 


27. 





be made from the coarse grid unless (perhaps) these values 
have been determined by a numerical method in which they are 
suitably affected by the fine mesh calculation. 

Specification 3 (exact solution) is a purely external 
specification, and is, therefore, computationally stable in 
the uniform grid. Figure 3 portrays this numerical result 
for the long wave case, and except for the "stair step" noise, 
the solution was qualitatively quite good. This is the first 
example of a specification which was stable in the uniform 
grid, but quite unstable in the meshed system. Computation- 
al noise initiated by slight inconsistencies between the 
numerical solution and specifications 3 - 8 in Table 1 pro- 


pagated upstream through the FMG in the meshed model. In 


AMPLITUDE 
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Fig. 3. Numerical solution (solid) compared with exact solu- 
tion (dashed) as t'= 10 for specification 3. See text and 
Table 1 for details. 
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the uniform grid system, this noise was reflected in a stable 
.manner by the inflow boundary, and does not cause exponential 
growth of these modes. However, in the meshed system, the 
computational modes represented an erroneous solution at the 
interface and thus contaminated the CMG solution, since the 
inflow of the FMG is the outflow for the CMG. The upper por- 
tion of Figure 4 depicts the behavior of the natural logar- 
ithm of system amplitude squared versus time during the in- 
tegration for specification 3 in both models. The solution 
in the meshed system remained stable until the computational 
modes reached the other interface boundary. At approximate- 
ly t'= 2,5 (see Figure 4), the amplitude began to grow ex- 
ponentially in the meshed model, while the solution in the 
uniform grid was accurate. 

Specification 4 (constant in time) was shown by Platzman 
(1954) to be computationally stable in a uniform grid, and 
this is verified in Table 1. The lower portion of Figure 4 
depicts the behavior of this solution for both uniform and 
meshed models. Because this specification is extremely in- 
consistent with the numerical solution, large amplitude 
computational modes were generated at the outflow boundary, 
and quickly filled the entire forecast region. The compu- 
tational modes grew linearly in the uniform grid model, and 
very quickly exponentially in the meshed system. Nitta 
(1962) mistakenly concluded that this boundary condition 
led to exponential growth in a uniform grid, but this was 


probably because his integrations were not carried out far 
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Natural logarithm of system amplitude squared vs. 
non-dimensional time for specification 3 (upper) 
and specification 4 (lower). See text and Table 1 
for details. 
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enough in time. One must be cautious in drawing conclusions 
from purely numerical experiments in which only a short in- 
tegration period is considered. Platzman (1954) has shown 
that exponential growth does not occur with such external 
specifications, and this is demonstrated by the linear growth 
depicted in Figure 4 for the uniform grid. 

Assigning the outflow point the value of the adjacent 
interior point (specification 3) excites rapid exponential 
growth of computational modes in both models, as portrayed 
in the upper portion of Figure 5, and indicated by the real 
part of the maximum eigenvalue in Table 1. It should perhaps 
be noted that each unstable specification ied in Table 1 
had many unstable modes, but only the maximum real part was 
listed in this report. The lower portion of Figure 5, how- 
ever, indicates that entirely different behavior took place 
in the uniform grid when the outflow point was assigned the 
value of the second interior point (specification 6). The 
modes amplified and damped quasi-periodically with no net 
growth, and this is indicated by the purely imaginary eigen- 
Values for this case. Computational stability for this 
Specification is a result of the centered spatial differenc- 
ame, in that the finite spatial difference at the first in- 
terior point was zero for all time. Spatial separation at 
alternating grid points is thus unable to affect the outflow 
boundary, as in the previous specification. The numerical 


solution, while computationally stable, is obviously useless. 
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non-dimensional time for specification 5 (upper) 


and specification 6 (lower). See text and Table 1 
for details. 
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Specification 7, on the other hand, in which the Tale at 
the outflow point is obtained by linear extrapolation from 
the two adjacent interior points, is strongly affected by the 
solution separation. Once initiated, rapid exponential 
Mouth occurs in both systems as can be seen in the upper 
panel of Figure 6. As noted earlier, one may be misled into 
concluding that the specification is computationally stable 
simply because the solution remains accurate for a short 
integration. As may be seen in Figure 6, the solution con- 
serves amplitude, and is quite accurate in the uniform grid 
for some time, but the onset of exponential growth is sudden, 
and the growth is explosive. 

The last case treated in these experiments was the util- 
ization of an upstream difference at the outflow point (spe- 
cification 8) similar to that suggested by Charney (1960). 
This is a dissipative differencing scheme (e.g., forward in 
time, backward in space), and in combination with centered 
space and time differencing at all other grid points, gave 
excellent results. This is also essentially the same speci- 
fication suggested by Shapiro and O'Brien (1970) for Hill's 
"telescope" model discussed earlier. Hill used an independ- 
ent' CMG forecast to obtain FMG boundary values (specifica- 
tion 2), and Shapiro and O'Brien showed that the use of up- 
stream differencing at outflow boundaries led to superior 
results. In the simultaneous integration of the meshed sys- 
tem, however, this specification also led to eventual expon- 
ential growth of the computational modes as seen in the 


lower portion of Figure 6. 
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DNA USE OF VARIABLE TIME INCREMENTS IN THE MESHED MODEL 

The stability analysis described in the previous section 
is valid only when the FMG is stationary and the same value 
of the time increment (At') is used in the numerical integra- 
tion of (4) in both grids. Use of a sufficiently small time 
increment to satisfy the Courant et al (1928) stability re- 
quirement in the fine grid not only reduces the accuracy of 
the solution in the CMG, but partly defeats one of the objec- 
tives of the meshed model, namely to reduce the computation 
time compared to a grid which is everywhere fine mesh. 

The next series of tests, therefore, were exactly the 
same as those discussed in the previous section, except that 
the time increments were in the same ratio as p, and higher 
wave numbers were used. Only those specifications which 
Pere computationally stable and produced accurate results 
were tested in this section. Thus, for the meshed model, 
only specification 1 (mutually determined) was used, and 
for the uniform grid model, only specifications 3 and 8 
(exact value and upstream differencing) were tested. The 
results for the uniform grid were, of course, as previously 
discussed, and the solutions for wave number one (not shown) 
closely resembled the upper portion of Figure 2. The meshed 
system was tested for p = 2 and p = 5, and also yielded a 
very accurate solution for both cases for wave number one. 
This was an extremely important result, as it verified the 


tacit assumption of linearity over the CMG time interval. 
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The FMG value at the external interface (j = 4 and j = 10 
in Figure l, for example) was determined after each FMG time 
step by linear temporal interpolation of the CMG solution 
for the tendency at this point. That is, the marching pro- 
cess required that the CMG time step be taken first, then 
the tendency at each FMG external boundary applied after 
each FMG time step. 

Higher wave number sinusoids were then imposed as the 
initial state to determine the behavior of the meshed grid 
when larger gradients passed through the stationary inter- 
face region. The first test was for an initial state of 
wave mabe: two. The meshed model numerical solution for 


p = 5 after 4000 time steps is compared to the exact solu- 


mom in Figure 7. Amplitude is predicted accurately, but 


AMPLITUDE 





Bor, Numerical solution (solid) compared with exact 
solution (dashed) for the meshed system with p = 5. The 
location of the FMG is denoted by the vertical dashed U 
For this test the FMG was stationary. 
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the phase error is approximately 0.03 wavelengths. The 
solutions for p = 5 were slightly more accurate than those 
for p = 2 (not shown) in every case, suggesting that higher 
values of p yield more accurate solutions over the entire 
grid. It will be shown later in nonlinear cases that the 
presence of the FMG in the model over strongly perturbed 
areas significantly improved the solution over the entire 
prediction domain, compared to a uniform coarse grid system. 
The uniform grid solution with the exact value specified at 
the outflow point (specification 3) for the wave number two 
case became extremely contaminated by computational noise, 
and hence it was eliminated from further tests. The upstream 
difference (specification 8), on the other hand, produced an 
extremely accurate کر ا نت‎ the uniform grid with vir- 
tually no error in either amplitude or phase for wave number 
two. 

The grid systems were then compared for the numerical 
solution of (4) when wave number four was imposed as an 
initial state. These solutions are compared in Figure 8, 
with the meshed model solution in the upper panel, and the 
uniform grid solution in the lower portion. The meshed grid 
solution in Figure 8 shows that the solution is accurate for 
the amplitude, but with a phase error of nearly 0.25 wave- 
lengths after 4000 time steps. This result may be compared 
to the control run, which was exactly the same grid without 
the FMG (cyclic continuity), where the phase error was 


slightly greater than that of the meshed system (not shown). 
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Fig. 8. Numerical (solid) compared with the exact solution 
(dashed). 


The upper panel depicts the comparison for the 

meshed system for wave number four with the location of the 
FMG denoted by the vertical dashed lines. The lower panel 
depicts the same comparison for a uniform coarse grid with 
Specification 8 as an outflow condition, 


and the exact value 
ase an inflow condition. 
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This result suggested that the presence of the FMG will not 
significantly improve the phase speed forecast of short wave- 
length features if they are allowed to pass through the inter- 
face. This is obviously a strong argument for a non-station- 
ary FMG which could remain centered over the short wavelength 
feature. The improvement of the phase speed solution in the 
meshed system over the control run was in approximately the 
same ratio as p. For example, if the FMG covered one-fifth 
of the total grid, then the phase error was approximately 
one-fifth less than the control run. 

The upstream difference (specification 8) yielded a rela- 
tively accurate solution, free of computational modes as may 
be seen in the tower panel of Figure 8. This solution cannot 
actually be compared to the one above it, because the exact 
value was specified at the inflow at each time step, but is 
shown mainly to demonstrate that this specification leads 
to good results, even with relatively short wavelength fea- 
tures. The result indicates that for uniform grids with 
centered space and time differencing, the upstream difference 
should be used at outflow points, This is in agreement with 
similar findings by Charney (1960) and Shapiro and O'Brien 
(1970). Further boundary condition tests with uniform grid, 
limited area models should be made, using nonlinear equations. 
One major disadvantage of the upstream specification is that 
System energy may not be conseryed. For very short range 
forecasts oyer limited areas, however, rather elaborate 


energy conserving boundary specifications and difference 
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schemes such as those proposed by Elsberry and Harrison (1971) 
may not be necessary. This was also suggested by Bushby 


and Timpson (1967) for limited area fine grid models, 
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IV. ONE-DIMENSIONAL MOVEMENT EXPERIMENTS 


A. MOVEMENT OF THE FINE GRID 

In general, the large scale flow has relatively small 
gradients of the dependent variables except in strongly 
perturbed regions. The results in the last section strongly 
indicated that to prevent large phase errors for short wave- 
length, large amplitude features, one should not allow such 
features to pass through the interface into the CMG. Hence, 
the FMG must be large enough to completely encompass the 
region of interest within the internal interface boundaries, 
and must be able to move in the appropriate direction if the 
feature approaches the interface. The movement tests showed 
that the meshed system satisfied two important criteria, 
namely that the solution remain computationally stable, and 
that it has comparable accuracy to that obtained on a uni- 
form fine grid. 

If the movement criteria, perhaps a region of large 
gradient approaching the interface was met, the FMG was 
moved one CMG spatial increment in the appropriate direction. 
This movement required new FMG values on the leading edge 
at both FMG time levels, and new CMG values on the trailing 
edge at both CMG time levels. It should be noted that the 
220.077600 ۶61605 time levels for the two grids are not tem- 
porally coincident. The trailing edge CMG values are simply 


assigned those of the previously spatially coincident FMG 
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points. Those at the current time level present no problen, 
but those for the previous CMG time level correspond to 
values from p previous time steps in the FMG. This time 
"alignment" is also a problem on the leading edge. On the 
leading edge, new values were obtained by first spatially 
interpolating the CMG values at both CMG time levels. The 
values were then interpolated in time to obtain previous 
time step FMG values, with the interpolation weighted by p. 
That is, the tendency at each point was only applied for 
one FMG time step. Thus, at the completion of the move, 
both CMG and FMG have values at appropriate points in time 
and space, and the integration continues, 

Another possibility of defining new values after the 
move was also examined. This involved the use of the values 
at the current time level only, which were obtained as des- 
cribed above. A Matsuno (Euler "backward") time step was 
then taken in the CMG and the tendencies determined the 
FMG external boundaries as before. After a similar restart 
in the FMG, the integration continued in the normal fashion. 
The only disadvantage to this method may be the undesir- 
ability of frequent restarts as would be the case in a 
rapidly moving system. Thus, any FMG movement discussed in 
the remainder of this report will be the former method of 


interpolation. 


B. LINEAR CASE 
In the previous sections, the mutually determined bound- 


ary conditions were shown to be computationally stable under 
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certain Conditions in the linear case. Examination of the 
use of different time increments in each grid while moving 
the FMG must be numerically treated. The non-dimensional 
form of (1) when scaled as previously described is 


9 9 
ee > ا‎ (5) 


In the numerical integration of (5), centered space and time 
0 ۱ ee 
differencing is used with ene 1 for each grid. “An initial 

state, which represents a sharp perturbation, or spike, 
within an otherwise undisturbed field, is given by the non- 


dimensional expression 


EO O CEOS 2 rele oie = (x'- x')] 


over one wavelength L with u = O elsewhere. This initial 
state is depicted in the lower left portion of Figure 9, as 
well as the analytic solution at non-dimensional times 
Meet and £' = 2, Accuracy of the numerical solution may 
be measured by comparing the amplitude syjuacred and phase 
velocity of the final numerical solution to the exact solu- 
Bion. 

As the time and space increments decrease, the numerical 
solution will tend to approach the exact solution. The 
upper left and lower right portions of Figure 9 portray 
numerical solutions on uniform fine and coarse mesh grids at 
t'= 2, and may be compared to the exact solution. The coarse 
grid (Ax =.01) produced an extremely distorted solution with 
a significant error in both phase velocity and amplitude 


squared. Trailing parasitic waves, or "wake," following the 


43 









SIME MESM REGION ág. 002 
COARSE MESH REGION AX. 01 


Ax+ 002 


AMPLITUOE A 
o 
o 


NON- DIMENSIONAL 
o 


oe EXACT 


Ax: 01 


۳ 


00 
o 1 2 
NON - OIMENSIONAL Time ) سے‎ 


Fig. 9. One-dimensional linear case. The exact solution 
for T = 2 is in the lower left panel. Numerical solutions 
for uniform fine grid, uniform coarse grid and meshed model 
are in the upper left, lower right and upper right respec- 
ively. Vertical dashed lines indicate initial and final 
positions of the FMG. 
spike disturbance are caused by truncation error, particu- 
larly with centered space differencing (see Wurtele 1961). 
The numerical solution is almost equal to the exact solution 
when the space increment is reduced by a factor of five, as 
seen in Figure 9, However, the "cost" of this improved 
solution is the utilization of five times the computer mem- 
ory and 25 times the computer time. 

If the two scales are combined as previously discussed 


to form a meshed system, a solution very much like that in 


the fine mesh uniform grid is obtained as may be seen in 
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Figure 9. During the course of the integration, 400 time 
steps were taken in the CMG, while 2000 were taken in the 
FMG. To maintain position centered over the spike required 
54 moves of the FMG toward the positive X direction. Ini- 
tial and final positions of the FMG are denoted by the ver- 
tical dashed lines in Figure 9. These linear experiments 
demonstrate an important aspect of the meshed system, namely 
the analytic treatment of stability given in the previous 
section may be extended to the case in which different time 
increments are used. Moving the FMG also had no apparent 


deleterious effects on the solution. 


C. NONLINEAR CASE 

The behavior of the meshed system in a non-linear envi- 
ronment is the next logical step in testing the technique. 
An interesting problem is the simple, but useful, example 
of hydraulic jump formation in a one-level hydrostatic fluid. 
This flow field has two characteristics which are particular- 
ly desirable for the test. First, Williams and Hori (1970) 
have shown that the solution is dependent on the space scale 
used, but in this case we can investigate the effect of re- 
ducing the grid increment only in the region of the jump. 
Second, the field develops sharp (eventually discontinuous) 
gradients in both pressure and velocity fields in a limited 
area. 


The non-dimensional equations describing the flow in a 


homogeneous inviscid fluid are, in flux form 
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tue to) d= -F Stu? +4)1-& +932 + alg tov (6) 


O 
2 Ivct +9) ]= -F luv (E +9)1-G +0)u (7) 
Lor o lu (54 9)) (8) 


Following Williams and Hori (1970), the equations have been 
L 


non-dimensionalized using a time scale T ہہ‎ , where 
وی‎ 2 
H is the mean depth of the fluid, the Froude number 
F = ےھ‎ , and the Rossby number R = —_ The hydraulic 
cn o fL 


jump forms due to the interaction of the non-linear terms in 


(6) and (8), and unless the Coriolis term dominates (i.e. 


ee 


R 
approximately t = 0.7. If the Coriolis term dominates, no 





>> F), the jump will occur at non-dimensional time of 


jump occurs and the equations describe an inertial gravity 
Mel lation in u and v. For thís test, both the Froude num- 
ber and the Rossby number were taken to be unity. Further 
details on an analytic treatment of the solution is given by 
Williams and Hori (1970), but will not be discussed further 


here. An appropriate non-dimensional initial state may be 


written 
u(x,0) = cos (21x) (9) 
$(x,0) = F [1/8 cos (4mx) + (52 - 1/8) cos(2mx)] (10) 
v(x,0) = 0 (11) 


and is depicted in the lower left portion of Figure 10. 
The numerical solutions of the energy conserving differ- 
ence equations (see Williams and Hori 1970) for fine, coarse 


and meshed grids are portrayed in the remaining sections of 
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Fig. 10. One-dimensional non-linear case. Initial state 
Mover left), and predicted fields of u, v, and ¢ for the 
indicated grid scales. The vertical dashed lines represent 
the initial and final locations of the FMG. 

Figure 10. As the gradient in the jump region increases, it 
becomes unresolvable in the coarse grid, and the solution 
quickly develops significant truncation error, resulting in 
a "wake" similar to that observed in the linear case. The 
fine grid, on the other hand, allows a superior numerical 
solution to be obtained, as noted in Figure 10. As in the 
linear case, the meshed model has virtually duplicated the 
solution of the uniform fine grid, although there is some , 
noise trailing the FMG. While remaining centered over the 


maximum gradient in this test, the FMG moved initially nine 


mes TG the left, then 17 times to the right. The initial 
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and final locations are again denoted by the vertical dashed 
lines in Figure 10. 

The one-dimensional tests of the meshing technique dem- 
onstrated several important points. First, computational 
stability based on the simple analytic case where the time 
increment was constant in both regions was demonstrated in 
more complex situations, where the time increment was differ- 
ent in each grid and the FMG was allowed to move. This is 
significant because use of time increments which are too 
small can not only reduce accuracy of the solution, but the 
Numerical integration is much slower. Second, the movement 
of the FMG did not appear to have any undesirable effects 
on the numerical solution in either the linear or the non- 
linear test. Movement of the FMG will be important for many 
potential applications, such as following tropical storms or 
other embedded smaller scale synoptic systems. Third, and 
probably most important, the meshed model produced a forecast 
quite similar to that produced by a grid which is everywhere 
fine mesh, while requiring much less computer time and 
memory space. As shall be seen in the next section, this 
is particularly significant in a multi-level model where 
many dependent variables are carried, and computer memory 
requirements may quickly rise to machine capacity if the 
grid is everywhere too small. Although the uniform coarse 
grid utilizes even less time and memory space, the forecasts 
in the region of large gradients are obscured by short wave- 


length noise. 
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The one-dimensional tests also strongly indicated that 
the mutually determined boundary specification in the meshed 
model is the only scheme which is computationally stable and 
produces usable results, at least for centered space and time 
differencing. Since it is generally desirable in geophysi- 
cal models to use such differencing, this is an important 
result. The use of an independently integrated CMG to spe- 
cify FMG boundary values was shown to be unsatisfactory, even 
for the long wave case. Use of such a specification in a 
complete non-linear model would certainly lead to aliasing, 
and possibly to non-linear computational instability. This 
Mas also pointed out by Phillips (1972). 

Another useful result of the boundary condition tests 
was obtained for the uniform grid. In this case, the up- 
stream differencing at the outflow point proved to yield 
superior solutions even for short wavelength, high amplitude 
features. Although the upstream differencing boundary con- 
dition may not conserve system energy, this aspect may be 
unnecessary for short term forecasts. In general circula- 
tion models of the atmosphere or ocean, however, non-conser- 
vation of energy may not be acceptable. This specification 
should also be tested in a non-linear model to more complete- 


ly determine the effect in such an environment. 
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V. TWO-DIMENSIONAL, NONLINEAR EXPERIMENTS 


An important step in the testing of the meshing technique 
is its application in a non-linear, multi-level prediction 
model. Since the ultimate objective of this research is the 
development of a meshed primitive equation model, results 
from a two-dimensional case will be useful for later exten- 
sion to the three-dimensional system. As in the one-dimen- 
Sional tests, simple analytic initial states will be imposed 
in the fine, coarse and meshed systems, and the forecasts 
from each will be compared. Boundary fluxes of heat and 
kinetic energy through the interface region will also be 
examined in both the meshed model and the corresponding por- 


tion of the fine grid model. 


Be THE PREDICTION MODEL 
The primitive equations in flux form governing zonally 


Symmetric flow in pressure coordinates are 





gu ے‎ = > E + fv (12) 


Is 


DP R/c 
٠ بو ؤ‎ ب٤‎ 


Q) 
I 


(16) 


where Q represents all sources and sinks of heat. Vertical 
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velocity in pressure coordinates (w) is carried at 100 mb 
intervals between 100 and 1000 mb in the ten-level model, 
while the remaining variables are carried at intermediate 


levels as depicted in Figure 11. Potential temperature and 
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Meee LL, Vertical distribution of dependent variables and 
pressure levels for the two-dimensional model. 
geopotential heights are also computed at 1000 mb for use 
in the hydrostatic and thermodynamic equations. The north- 
south boundaries are insulated free slip walls, while the 


lower boundary condition is 


d¢ 
1000 
0ج2‎ = Q (17) 


Vertical velocity is obtained by integrating (14) downward 
from 100 mb where w is assumed to be zero. The time inte- 
tion is centered after an initial forward time step, with 
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certain differences in the meshed system which will be dis- 
cussed later. Other aspects of the basic numerical model 
are given by Elsberry and Harrison (1971) and will not be 


discussed further here, 


B. TWO-DIMENSIONAL TESTS 

An initial state representing a zonal trade wind regime 
of low-level easterlies replaced by westerlies aloft was 
achieved by imposing an appropriate horizontal temperature 
gradient. This initially undisturbed pattern, depicted in 
the vertical cross section in Figure 12, is distorted as a 
constant, analytic heating distribution is applied, within 
a relatively narrow band. The heating function was similar 


momthat used by Holton (1971) in diagnostic experiments, 





PoP E 
yey, A ap.) TIPS 
Q = Q exp [- ( o) e | Po Pmax P-P nax maz (18) 
P 
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Here 7 and y are the locations in the vertical and hori- 
zontal where the heating function is a maximum. The para- 
meter a determines the relative latitudinal extent of the 
heating pattern and B determines the rate at which the func” 
tion decreases above the maximum vertical level. The verti- 
cal section of the heating plotted in Figure 12 results when 
the following parameters are chosen 
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ES“ 12. Initial wind field (solid) in E pressure. 

The heating pattern (dashed) is a maximum at 550 mb, with 
Molines labeled in K dayl. 

This heating pattern simulates the heat released by a zonal- 
ly symmetric band of convective elements. To offset this 
heating, radiative cooling over the region was chosen such 
that the net heating at each pressure level was zero. In 
the coarse grid, most of this heating will be applied at the 
point y = A (see Figure 12), while in the fine grid with 
one-fifth the space increment, the heating will be smoothly 


applied over many grid points. The sum (Q) of the imposed 
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heating and compensatory radiational cooling in (14) is 
chosen so that the volume integral of heat energy is constant 
Mm time. 

Figure 13 depicts the 24 hour forecast field of 1000 mb 


heights for each of the three models. For this forecast, 
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meee ts. Predicted heights at 1000 mb (solid) for the 
models indicated. 480 time steps have been taken in the 

CMG and uniform coarse grid, while 2400 have been taken 
mthe FMG and uniform fine grid. 

480 three minute time steps were taken in the coarse mesh 
region (AX = 45nm), while 2400 time steps were required in 
the fine mesh regions of the uniform and meshed models, 
because of the factor of five in the spatial increment. The 


effect of the localized heat source in the uniform coarse 


grid is quite evident as can be seen in the lower portion 
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ET Ticure 13. Because nearly all the thermal forcing is 
applied at one grid point (column) in the coarse grid, a 

2Ax mode is excited. The corresponding forecasts by the 
meshed and fine grid models are shown in the center and 
upper portion of Figure 13 respectively. In the area of 
maximum heat release of both systems, nearly identical 1000 
mb height fields are predicted. Even in the CMG of the 
meshed اٹ‎ the forecast is quite similar to that plotted 
"oT every fifth point within the uniform fine grid system. 
Perhaps a more sensitive indicator of the performance of the 
meshed model is the response of the vertical motion field to 


the imposed heating pattern. Figure 14 depicts this field 
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LAN ertical section of the 24 hour vertical motion in 
mb sec7lx104 for the FMG (solid) portion of the meshed model, 


and the corresponding region of the uniform fine grid model 
(dashed). 
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morschesneshedrand uniform fine grid within the FMG region 

of the meshed model. Very little vertical motion occurs out- 
side this limited area and is therefore not portrayed. A 
deep region of ascent is produced by the heating, with weaker 
compensating descent on either side. In the meshed grid, the 
ascent cell was stronger and slightly wider after 2400 FMG 
time steps, probably due to non-linear effects at the inter- 
face. The heat gains within the FMG meshed region must be 
transferred horizontally because of the requirement for heat 
balance at each level, and boundary fluxes across the inter- 
face are an important consideration. 

The accumulated boundary fluxes of kinetic and thermal 
energy through the interface in the meshed system, and the 
corresponding area in the uniform fine grid, may be obtained 
by integrating these fluxes in the vertical and around the 
entire area after each FMG time step. In Table 2, this com- 
parison is made for kinetic energy between the two grid 
systems. The row labeled "inner" refers to the FMG portion 
of the meshed model, and the corresponding area of the uni- 
form fine grid system, with the remainder of the grid area 
represented by the "outer" row of data. The larger genera- 
tion in the meshed system is due to the somewhat stronger 
ascent cell, but despite this, the boundary flux of kinetic 
energy was the same in both models. The dissipation quantity 
for the meshed model is a result of a horizontal momentum 
diffusion term used in this system and will be discussed 


rer To dissipative terms are included in the other 
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Table 2. Kinetic energy balance for 
meshed and uniform fine grid models 
(joules cm 2 day Lys. 


ATF E Generation | Dissipation | Boundary Flux 


Inner +2.08 +2.54 -0.46 


Fine Grid 


Outer +1.51 +1.05 


Inner +2.38 EI -0.93 





Meshed Grid 


Outer +1.25 +0.79 


+0.46 


formulations, thus the balance is mainly between boundary 
flux and generation. The qualitative comparison of the 
kinetic energy balance is good, with the total kinetic ener- 
gy change over the entire forecast region almost identical 
in both systems. 

Because the radiational cooling exactly offsets the 
imposed heating, the net change in thermal energy in Table 3 
is necessarily zero over the entire forecast region. The 
slightly larger boundary flux of thermal energy in the 
meshed system is again due to slight inconsistencies at the 
interface. As in the case of the kinetic energy balance, 
the qualitative comparison of the models in Table 3 is excel- 
lent, with less than three percent difference in the boundary 


fluxes. 
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Table 3. Thermal balance for meshed and 
uniform fine grid models (deg.pot.temp.cm day 
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C. POMEPUTATIONAL DEVICES USED IN THE MULTI-LEVEL MODEL 














Cooling {Boundary Flux 


Fine Grid 






Meshed Model 





Early tests not discussed in this paper were conducted 
to examine the behavior of the meshing technique in the 
primitive equation model. In these tests, two aspects re- 
lated to the interface boundary conditions and the centered 
Space and time differencing became evident. Spatial separa- 
tion in the FMG resulting from slight inconsistencies in 
the imposed boundary condition created the first problem. 
Although the solution on the FMG near the interior coinci- 
dent point is "nearly exact," the CMG solution at the exter- 
nal interface is not as accurate, and will differ somewhat 
from one obtained on a uniform fine grid, even though the 


integrations. are simultaneous. The tacit assumption of 
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linearity over each CMG time step also adds to the incon- 
sistency. Once spatial separation existed within the FMG, 
setting the CMG boundary values over an odd number of steps 
triggered high speed internal gravity waves which propagated 
from the interface into the CMG. When reflected from the 
north-south walls, the gravity waves tended to interact with 
those propagating in the opposite direction and growth on a 
small scale was observed, eventually leading to non-linear 
computational instability. The minor spatial variations at 
the interface were eliminated by applying an area (line) 
average of the surrounding five points, rather than the 
point value at the single coincident point. As a final con- 
straint on the spatial separation in other regions of the 
FMG, a horizontal momentum diffusion term of the Laplacian 
type was employed in this portion of the model. A relative- 
ly small coefficient of 10° ae ٦ was sufficient to control 
the separation. 

Even with the spatial separation controls discussed 
above, the slight inconsistencies in the boundary values 
derived from the CMG solution caused the solution in the FMG 
to begin to separate in time. This was a small effect, but 
after several hundred steps in the FMG, it was evident that 
control was needed. Thus, Matsuno (Euler "backward") re- 
starts at experimentally determined intervals were periodi- 
ails applied in the FMG and the time separation was elimi- 
nated. These controls on the meshed system appear to have 


very little effect on the accuracy of the solution as 


59 





evidenced in the last section. While it is true that higher 
order differencing and/or interpolation schemes might elimi- 
nate the need for these simple controls, the additional com- 
putation time required for such formulations would probably 
negate their value. Even with the centered space and time 
differencing, the time increase for the integration in a 
three-dimensional meshed model will not be trivial. 

It is clear from examination of Tables 2 and 3 that the 
meshing technique conserves total thermal and kinetic energy 
when the finite difference equations are written similarly 
to the energy conserving flux form described by Arakawa 
(1966). This is important not only because of computational 
meabatlity Criteria, but also for long term application. 
General circulation models of the atmosphere and ocean often 
require increased resolution in certain areas, and it is 
crucial that no spurious sources or sinks of energy are 


introduced by modifications to the grid network. 
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e REEE DIMENS TONAL EXPERIMENTS WITH A 
ERIMITIVE EQUATION MODEL 

The one- and two-dimensional tests demonstrated the com- 
putational stability of the meshing technique in a variety 
of environments. Nonlinear two-dimensional cases illustrated 
that realistic energy exchanges occurred at the interface, 
and that improved solutions were obtained with the meshed 
grid system. The three-dimensional cases discussed below 
will demonstrate that these features of the meshing technique 
are maintained in more complex situations. 

Important aspects of three-dimensional meshed grids are 
initialization of the fine mesh region, physical size of the 
fine mesh region and ratio of grid increment reduction. The 
latter two considerations are strongly interdependent for 
three-dimensional models, since from computer memory and time 
viewpoint, increasing the physical size of the FMG restricts 
meerd scale reduction. On the other hand, the FMG must 
be large enough to completely encompass the region of inter- 
est, but still have a ratio of reduction to a small enough 
spatial increment so that the FMG can allow development on 
the scale desired. Thus, one must evaluate the type and 
scale of the feature of interest, and design the FMG with 
mese IN mind. For example, the same FMG specifically de- 
wened for tropical storm studies would probably be unsuit- 


acle for a frontal scale feature. In the latter case, an 
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FMG with larger physical size, but with a smaller scale re- 
duction would be more appropriate. 

An extension of the concept of optimizing FMG design to 
suit the scale of the feature of interest is to utilize a 
EraGdual reduction in grid scale. This would essentially 
allow one to "focus" on a region of interest with a very 
fine grid, while still maintaining relatively fine resolution 
in the surrounding area, Certain difficulties and limita- 
tions associated with this type of model arise from applica- 
tion of the mutually determined boundary specification, and 
will be discussed in detail in a later section. A primary 
advantage of the gradual reduction model is the reduction 
of computer memory and time required to obtain improved solu- 
tions over the same area as an abrupt reduction to the 
smallest scale. 

As in the earlier tests, simple analytic initial states 
were used in the three-dimensional models to idealize the 
initial conditions. This idealization is desirable since 
the primary objective of these cases was to demonstrate the 
meshing technique in the single reduction (SR) and multiple 
reduction (MR) models, observe the effect on the large scale 
solution, and compare the results to those obtained on a 
uniform coarse grid. Improvements to all the versions of 
the model would have to be accomplished through real data 
forecasts, with subsequent statistical verification. There 
are, for example, many possible ways to initialize the FMG, 


and many real data tests would need to be made to determine 
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the optimum method, Other facets of the meshed models, such 
as diffusion coefficient magnitude, frequency of restarts, 
etc., have only been cursorily examined in this report and 
could probably be optimized through further numerical test- 


ing. 


CE THE PREDICTION MODEL 

The basic model utilized in all three-dimensional exper- 
iments was similar to that described by Harrison (1969) and 
Elsberry and Harrison (1971, 1972), except that fewer levels 
in the vertical were used. As in the two dimensional case, 


the equations written in flux form are 


2u =- L(u) + fv - m 2È (19) 
ay = - L(v) - fu - n $$ (20) 
<2 = - 100) +Q ao 
1۰ب‎ + (22) 
a دلاو‎ = 


where Q represents all sources and sinks of heat and 


2.0 su ð Sv d 
es a + a ا‎ ate a (wS) 


L(S)= m 
ior a general variable S. The map factor is m, with the grid 
SS Mercator projection true at 22.5N. The vertical distri- 
bution of variables and pressure levels in the three-dimen- 
sional version is depicted in Figure 15, with the upper and 


lower boundary conditions identical to those in the two-: 


dimensional model. The north-south boundaries of the three- 
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lS. Vertical distribution of dependent variables and 
pressure levels for the three-dimensional models. 
dimensional model are insulated, free-slip walls, while 
eyclic continuity is assumed in the east-west direction. 
As in the two-dimensional case, the difference equations 
were written similarly to the energy conserving form des- 
cribed by Arakawa (1966). 

An example of a horizontal grid level for the SR meshed 
system is depicted in Figure 16. The CMG indices shown are 
acronyms by which the location of the FMG may be determined. 
For example, IWE represents index I(east-west), West side, 
and External interface. Thus, if IWE and JSE are known, the 
location of the FMG is uniquely determined, and these, of 


course, are time dependent as the FMG moves. To avoid 
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ES l6. An example of a horizontal pressure level for the 
three-dimensional meshed model. For the example shown, the 
TIC increment is one-fourth that of the CMG (p = 4), and 
E TIC covers Six CMG intervals in each horizontal direction. 
interaction with the CMG boundaries, the FMG is not allowed 
Eo move within three CMG grid intervals of the walls, or 
the east-west extremes of the grid. For these tests, the 
CMG spatial increment was 180 nm, and for the SR model, the 
FMG increment was 36 nm (p = 5). Because of excessive com- 
puter memory and time requirements, no uniform fine grid 
comparisons were possible for any of the three-dimensional 


experiments. Other than the restrictions above, the location 


PELE EMG Shown in Figure 16 is quite arbitrary. Details 
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of the MR model will not be discussed in this section, and 
ehe first series of tests to be described utilized the SR 
model with p = 5. 

The CMG used for all three-dimensional tests was 25 x 20 
in the east-west and north-south directions respectively, as 
may be seen in Figure 16. In the SR meshed model, the FMG 
covered 6 x 6 CMG grid intervals, or a 31 x 31 FMG. Because 
p = 5 in this model, the time step in the FMG was one-fifth 
that in the CMG, and it is quite evident that the FMG calcu- 
lations required most of the computer time, since the number 
Of total calculations increases as the cube of p. This again 
ted the desirability of keeping the FMG just large 
enough to encompass the region of interest. It was experi- 
mentally determined that each time step took approximately 
eight times longer in the meshed system, as opposed to the 
uniform coarse grid. This may be compared to the time in- 
crease of about 125 times for a uniform fine grid over the 
entire region of the coarse grid model. 

Initial location of the FMG is achieved by identifying 
the desired location and appropriately defining IWE and JSE 
(see Figure 16). At the interior Boe ene points (ISI, 
JSI, for example) the CMG boundary values were determined 
by an area average of the surrounding 25 FMG values before 
۳۰٠۰٠٠٠٦٠٣٦٦٦٦6 >٠ ٥م‎ ,- For convenience in the actual coding, 


PementceLmteCrlOn coincident points were replaced by the area 





average, realizing that no forecasts are made at the CMG 


points which lie within the FMG area. After each complete 
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time step in the CMG, the tendencies around the external in- 
terface were spatially and temporally interpolated. As in 
the earlier experiments, the FMG external boundaries were 
then updated after each FMG time step. Movement of the FMG 
was accomplished in exactly the same manner as in the simp- 


ler models. 


DEN TTIALIZATION OF THE MESHED GRID SYSTEM 

Initialization of the FMG is an important point to consi- 
der in the meshed model. Since mesoscale data are seldom 
available, there are two basic alternatives, that is to 
either attempt a fine grid analysis or to interpolate data 
from the CMG. Hill (1967) has compared fine grid analysis 
interpolation for initialization of the FMG, and noted 
that if analyzed data are used, the effects of errors in 
height or temperature gradients increases as the grid dis- 
tance decreases. It would seem logical, therefore, that 
interpolation from the CMG is the only real alternative, un- 
less perhaps an objective analysis scheme was used, with 
basically the same order interpolation as the interpolation 
from the CMG. If mesoscale data are available, however, one 
might consider modifying the interpolated values. For 
example, if the interpolated data were to undergo a diagnos- 
tic phase such as the balance equation, the observed meso- 
scale data could be incorporated before the diagnostic phase, 
thereby enhancing the representativeness of the analysis. 
Since no small (sub-grid) scale data can be added by inter- 


polation, one must hope that the elements of the small scale 
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will naturally develop during the integration as a result of 
the large scale forcing or parameterized fluxes of heat, 
moisture and momentum (e.g. Bushby and Timpson 1967). This 
is a fundamentally crucial concept to fine mesh forecasting, 
for if it is not true, grid increments need be no finer than 
the distance between observation stations. 

After the FMG data has been interpolated from the CMG, 
several alternatives exist. That is, the values of the 
dependent variables may be used as is, or undergo some type 
of diagnostic phase. Wang and Halpern (1970) have compared 
the use of pure interpolated data versus that which had been 
balanced by means of the linear balance equation, and obtained 
better results with the interpolated data. Their model was 
barotropic, and smoothers were used, so it is datada to 
evaluate their results as applied to the multi-level primi- 
tive equation model. Thus several experiments were conducted 
us evaluate the optimum initialization for the analytic 
initial states. The first attempt was to merely derive ali 
values of the dependent variables from the CMG to the FMG 
by lineal interpolation, and start the integration with no 
diagnostic phase, Almost immediately, large amplitude very 
short wavelength gravity waves were excited, and aliasing 
very quickly led to non-linear computational instability. 
Pure interpolation was clearly unsatisfactory for the meshed 
model, at least when only small diffusion was employed, and 
there was no implicit mechanism for selectively damping high 


frequency modes. 
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A consistent diagnostic equation which uses winds to 
rive heights may be obtained by forming a divergence equa- 
tion from (19) and (20). Such an equation is particularly 
applicable where the Rossby number U/fL is not small, such 
as in the tropics, or for small scale disturbances. The 


equation may be written 


Ru dv 0 0 IV du 
o ay) ox PO gy Lv) +E T ay 
of om 0 2 
- u 3 - E - mV Q (24) 


This divergence equation as shown is not exactly consistent 
with the form of the continuity equation (22), but because 
the map factor was taken to be unity for these experiments, 
the additional multiplications by m are not shown. In the 
first SR model test to be described (tropical storm simula- 
tion), the winds at all levels were analytically imposed. 
The divergence equation (24) was then solved by sequential 
Byerrelaxation to obtain the geopotential heights, and hence 
temperature, with the assumption that the time change in 
divergence is zero over one time step. For this case, the 
original winds were non-divergent, but it should be noted 
that this is not necessary for real data or interpolated 
data. In the second SR test (frontal simulation), the FMG 
was initialized quite differently. In this case, data were 
obtained from a forecast on the uniform coarse grid model. 
nS dala was thus the result of several model days of inte- 
eration, and hence the vertical motion fields were well 


io Dedo Winds were interpolated from the CMG to the FMG, 
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and a consistent vertical motion field was obtained from (22) 
and so the nonlinear terms in (24) were complete. This 
approach was the last one tested during this research, but 
seems to show promise as a possible technique when real data 
is used to initialize the model. In all SR and MR experi- 
ments, the CMG was initialized with the use of (24), and 
boundary conditions for this grid are those suggested by 
Elsberry and Harrison (1971) and will not be discussed here. 
SR TEST 1: A SIMULATED TROPICAL STORM WITH 
DIABATIC HEATING 

An important meteorological application of nested grid 
Systems is in tropical prediction models. Formation, devel- 
opment and subsequent movement of tropical storms and 
Similar disturbances could be more accurately predicted by 
a properly located FMG. The first three-dimensional test 
to be discussed involved the SR model, and an analytic simu- 
lation of a tropical storm scale disturbance. As noted 
earlier, analytic initial states were used throughout in 
wder to idealize the initial conditions for the tests of 
the meshing techniques. Forecasts of the movement and devel- 
opment of the simulated disturbance were then compared for 
the meshed and uniform coarse grid models. Boundary fluxes, 
enee nergy balance and temperature patterns are treated 
Is testo as well as several other aspects of 
Ene integration. 

The storm simulation was kept as simple as possible to 


isolate possible effects of the meshed system on the 


70 





numerical solution. An isolated low pressure center in the 
trough of a wave in the easterlies may be simulated by speci- 


fying the winds uniformly at all levels by 








90 
1 0 
U = - ہے‎ — (25a) 
f oy 
96 
urn 2- —_ (25b) 
0 
where 
x-xX 
YE = F(y)+ H, sin(k,x)sin (k,y)- H, exp ( = ©) exp 
A 
( B ) (250) 


Since the winds are constant in the vertical, the initial 
disturbance is a circularly symmetric, non-divergent vortex, 
The expression for 9, is convenient, because the size, shape 
and location of the resulting balanced low center are easily 
controlled. The independent variables in (25) have been 
normalized by the length and width of the domain for simpli- 
city. The magnitude of the basic zonal flow was approximate- 
IS 4.5 m ae and is determined by F(y), a linear function 
of latitude. The amplitude of the wave in the easterly flow 
is governed by Hos while the central value of PL is deter- 
mined by Hz. The parameters a and B govern the longitudinal 
and latitudinal extent of the depression center, and ky and 
ky are the domain wavelengths in the east-west and north- 
south directions respectively. The low center may be appro- 


priately located in the trough by specifying x, and Yo to 


0 
Coincide with the minimum values determined by the middle 


Bern Jin (250). 
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Be. 17. ۳5۰٣٠٠٠٦٦٠٠٦6 ا‎ field (m) 1400m at 850 mb after 
Meer alization procedure described in text. The initial 
tion of the FMG is denoted by the dotted outline. The 
Mana uniform coarse mesh grid scale is indicated in the 
lower left corner, 


eure 17 depicts the resulting balanced height pattern in 


MCMC for the following choice of parameters: 


H, = 150 2 a 
O a 
3 

a = 1 

B = 1 
= 

Ep = 08 
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After the CMG balanced geopotential heights were obtained, 
each of these parameters, as well as L(u) and L(v) in (24) 
were interpolated around the external interface. These 
values were held fixed during the subsequent diagnostic phase 
of the FMG, and thus assured a consistent merging of the two 
eres in the interface region of the FMG. Specifying the 
geostrophic wind analytically in the FMG creates a simulation 
of enhanced data in this region. The values at interior co- 
incident CMG grid points were exactly those which would have 
been obtained by interpolation, but there was no assumption 
of linearity between these points. In an extremely idealized 
sense, this had the same effect as subjective modification 

of the initially interpolated FMG data with mesoscale data 
obtained from some other source. 

ET portion of the 850 mb geopotential height field 

within the dotted outline of the FMG in Figure 17, may be 
compared with Figure 18, eon depicts the same field result- 
ing in the FMG when the above initialization is carried out. 
Near the interface there is a smooth blend, because of the 
boundary conditions discussed in the preceding paragraph. 
In the interior of the FMG, however, there is more detail in 
The initial state. Addition of such data would clearly im- 
prove the forecast in the FMG, and perhaps simulate develop- 
ment of features on this scale. 

eee Geom Che Mropical’ Storm Simulation (SR Model) 

Development of the simulated storm was artificially 


imduced in a manner similar to that used in the two- 
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Fig. 18. Balanced height field (m) -1400m at 850 mb in the 
Ib tollowing the initialization described in the text. The 
86ا۶‎ erid scaie is indicated in the lower left corner. 


dimensional test. In actual tropical storms, one expects to 
find maximum latent heat release in approximately the middle 
portion of the troposphere, thus an analytic heating function 
was used to simulate such heating. This function had the 
advantage that it could easily be positioned in both the 


vertical and horizontal and is given by 








xx N 
Q = Q, expl- e] exp ٭ رڈرقتے۔‎ (26) 
Pp TP > 
ee Ten P < Pmax 
L p =P 
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where xo» Yo and Prax are ٣٢۷٣٠٠٦٦٦272 0٦ ۲۹1۸ ۹718 vertical locations 
of the heating maximum. The parameters a and B determine the 
longitudinal and latitudinal extent of the heating. The rela- 
tive lack of vertical resolution in the model caused more than 
half of the heating to be applied at Bee To simulate heat- 
at levels between 550 mb and 250 mb, the heating at 250 mb 
was set equal to that at 850 mb. As in the two-dimensional 
case, a radiation profile was derived such that the net. heat- 
ing at each pressure level was zero for the entire spatial 
domain. 

The heating pattern for the FMG at ےت‎ Ls depicted 
in Figure 19 for the following parameters: 


Qo = 5 X noes K 5 


Qa = 1 

B 21 
ج7‎ 550 mb 
Po = 1000 mb 


Meee ne distribution in the uniform coarse grid 

(not shown) was slightly different because the area averages 
of the FMG heating were applied at coincident points, but 

the net heating over the entire spatial domain was the same 
Tor both models, In both uniform coarse grid and SR models, 
the horizontal heating maximum was linked to the location 

Bi The 1000 mb Low center. In the FMG of the SR model, the 
location shifted with the 1000 mb low center in intervals of 
one FMG spatial increment. Thus in the FMG, the heat source 


was always in approximately the proper location to enhance 
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feo 19, Heating distribution (K کہ‎ in the FMG at Da 
DD). The distribution at the other levels was 
similar, but with smaller values. The FMG grid scale is 
indicated in the lower left corner. 

Bie generation of kinetic energy. atenas manner, the link 
between boundary layer ascent and latent heat release (con- 
onal instability of the second kind, or "CISK") was 
roughly simulated. When the FMG moved (one CMG increment), 
the heating maximum was redefined over the new location of 


the low center in the FMG. This small scale detail was, of 


Bourse, not possible in the uniform grid model, where the 
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heating maximum location was shifted in increments of one 
CMG spatial interval, again dependent on the 1000 mb low 
center. 

It is apparent that proper location of heating with 
respect to the depression center will improve the forecasts 
of both development and movement, at least in a fine grid. 
This may not be true, however, in uniform coarse grid models. 
Elsberry and Harrison (1972) showed that the use of differ- 
ent parametric latent heating schemes had little effect on 
the predicted movement of hurricane Betsy (1965) in a uni- 
form coarse grid, ten-level primitive equation model. Al- 
Eroueh Elsberry and Harrison's model was very complete in 
that friction, latent and sensible heating and radiation 
were included, the movement can only be subjectively deter- 
mined by "guessing" at the location of the low center at 
the completion of the forecast. Thus as the grid scale 
decreases, the location may be more accurately determined, 
thereby yielding a more valuable movement forecast. Also, 
as the scale decreases, the parametric relationships pro- 
bably become more valid, since the underlying assumptions 
have greater physical 08 on a smaller scale. It seems 
reasonable to expect higher verification score in real data 
cases, when the same parametric schemes are employed in the 
meshed system, as opposed to a uniform coarse grid model 


such as that used by Elsberry and Harrison (1972). 
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2 EES Oi ChE ETE Tropical Storm Simulation Case 


Two significant results of the tropical storm simu- 
lation in the SR model involved the movement and development 
of the low center. The phase error of short wave systems is 
strongly linked to the finite grid scale, and, as seen in 
earlier linear experiments, the meshed model reduces the 
error. Reduction of phase error is particularly crucial in 
we dicting the movement of tropical storms, where accuracy 
is extremely important. 

Te piediCcted 46 hour, 8530 mb height field for the 
۲٣۰٠۰ orm coarse grid is depicted in the upper portion of 
Figure 20, with the same forecast for the meshed model shown 
in the lower panel. The meshed model result was obtained by 
Ne the CMG prediction, but replacing interior coin- 
cident values by those from the FMG forecast. Since the 
۰٠٠۹٦ forecast location of the low center in the FMG was 
between CMG grid points, this location was, therefore, taken 
from the FMG forecast. As may be seen in Figure 20, the 
meshed model predicts greater westward movement of the low 
center than the uniform coarse grid model. Since the mean 
zonal flow velocity was approximately 4.5 m ae a pure 
advective westward movement of about 775 km would be expec- 
medwauring the 48 hour forecast. Location of the low center 
in the uniform coarse grid is rather subjective, and thus 
movement is difficult. to ascertain with accuracy finer than 


the grid scale itself, but Table 4 summarizes the two main 
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Peewee cedtiction of 350mb heights (m)-1400m at 48 hours. 
ep per portion is the uniform coarse grid forecast, the 
lower that of the meshed model. See text for details. 


79 





Table 4 


Comparison of the uniform coarse grid and 
meshed model 48 hr height forecasts at 
850 mb (low center) 





Model Movement Development 


Distance Forecast | Change 
(km) (m) (m) (m): 
| 


facets of the prediction, that is movement and development. 


Uniform 
coarse 





The meshed model presented distinct advantages in evaluating 
the forecast fields of the dependent variables. An example 
Si the detail available is portrayed in Figure 21, the 48 
hour forecast 850 mb heights. Location of the low center 
was accurate to the scale of the FMG, thus both location and 
direction of movement were more accurate. This information 
was then transferred to the large scale analysis shown in 
the Lower portion of Figure 20, which is perhaps more con- 
venient for display purposes. Thus the additional detail 
available from the FMG forecast is incorporated into a large 
area analysis in a consistent manner. The initial and final 
Bo itionssof the FMG are denoted by the dotted outlines in 
Figure 17 and Figure 20, respectively. The FMG moved toward 


eC vest twice during the 43 hour forecast, and as seen in 
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Mme D1. Prediction of 850 mb heights (m) -1400m in the 
fine mesh after 48 hours. The FMG has moved twice during 
the forecast period. The FMG grid scale is indicated in 
the lower left corner. 
Figure 21, the low center was again approaching the western 
side of the FMG. The large scale flow patterns were pre- 
تہ‎ ۰ ۹3 quite Similarly in both models. 

Another example of the detail available in the FMG 
forecast is shown in Figure 22, the 48 hour vertical (pres- 


ure) velocity pattern at 700 mb in the FMG. The low center 


at 850 mb is denoted by the bold C, and areas of maximum 
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meee 22. Vertical (pressure) velocity (10 “mb sony in the 
Micearter 48 hours at 700 mb. The location of the 850 mb 
low center is denoted by the bold C, and the FMG grid scale 
is indicated in the lower left corner. 

ascent are shaded. Although the heating is symmetric with 
respect to the low center (see Figure 19), there is a ten- 
dency for spiraling of the narrow bands of maximum ascent. 
Ne spatial inecrement of 36mm in the FMG prevents any very 


small scale details of the tropical storm simulation to be 


Motera edo but the overall pattern in Figure 22 is roughly 
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suggestive of observed data. If the model employed a para- 
metric latent heating scheme dependent on vertical velocity, 
the shaded regions of ascent in Figure 22 might then be in- 
terpreted as cloud bands. The fact that such regions of 
ascent form in the model, without friction or a sophistica- 
ted latent heating scheme, suggests that such patterns may 
be independent of friction, and are probably linked to the 
dynamics of the flow itself. 

Kinetic energy balance was also calculated for the 
CMG and FMG of the meshed model, and for the uniform coarse 
grid model, and is summarized in Table 5. As in the two- 
mensional case, horizontal diffusion of momentum was em- 
ployed in the FMG of the meshed model to suppress short 
wavelength computational noise, and therefore this region is 
the only one with dissipation included as a term in the 
balance, There are three separate balances because movement 
of the FMG required that a new balance be started at the 
completion of each move. Early in the prediction period 
Tle 54) the diabatic heating had not yet caused the 
development of strong ascent regions in the FMG, and large 
dissipation governed the reduction of kinetic energy over 
this period in the meshed model. Generation of kinetic 
energy over this period was even less in the uniform coarse 
grid model, but because there was no dissipative mechanism 
ar Boundary flux, all the generation went to increase the 


total kinetic energy. 
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Table 5 


a. Kinetic Energy balance for an and uniform coarse 
grid models (10—“joules cm *day~') up to first move CMG 
time step 0-80. 


Boundary 
Model AKE Generation Dissipation flux 


Meshed 
CMG +3.2 +1.5 +1.7 
Uniform 
grid +4.5 +4.5 


b. Kinetic energy balance for meshed and uniform coarse 
grid models up to second move CMG time step 80-160. 





A Boundar 
Generation DpisSipation A 2 
Tr us 


A A کچ‎ 
ne 


Meshed 
er 


Uniform 
grid 





CG’ Kinetic energy balance for meshed and uniform coarse 
grid models up to 48 hr CMG time step 160-240. 


Boundary 








Model AKE Generation Dissipation 
flux 
PNG as | 13.5 ESA 
Meshed 
CMG 
Uniform 
grid 000 7)0 - | ae 





In the intermediate period (Table 5b), the meshed 
system nearly restored the kinetic energy of the meshed sys- 
tem to the initial value, and the net generation was larger 
than the dissipation. Development of the regions of strong 
ascent as depicted in Figure 22 were established during this 
period, and the boundary flux across the interface remained 
approximately constant. Very little change occurred in the 
wnt orm grid model in this phase of the prediction period. 

During the last phase of the prediction (Table 5c), the 
vertical motion pattern depicted in Figure 22 was fully 
established, and active over the entire phase in the meshed 
model, with accompanying large generation of kinetic energy 
in the FMG. The two-dimensional experiments previously 
indicated that the mutually determined boundary conditions 
wened efficient interface transfer of kinetic energy, and 
this result was extended to the three-dimensional model as 
well. This may be seen 1n Table 5c, in that the greatly 
increased generation in the FMG did not cause excessive 
increase in the FMG kinetic energy, instead, FMG kinetic 
energy increased only slightiy over the initial value, and 
the boundary flux easca with the increase in generation. 
The CMG apparently converted the kinetic energy to potential 
energy by gradually subsiding the warm upper level outflow 
from the interface over a large area. 

Thermal balance in the two-dimensional test was ob- 
tained by performing all model computations in "double 


Poeeetciom On the IBM 360/67 at the Naval Postgraduate 
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School. This was necessary because the normal seven signifi- 
cant decimal places were not sufficient for balance, due to 
the magnitude of certain sums and forcing functions in the 
predictive equation (21). Increasing the decimal signifi- 
cance by a factor of two, however, nearly doubles the computer 
memory space requirements, and thus was impossible for the 
three-dimensional model. Behavior of accumulated heating and 
thermal boundary fluxes in time were examined, however, and 


these comparisons are portrayed in Figure 23. The accumulated 


Accum Flux 
Accum Hig 





Fig. 23. Comparison of accumulated thermal boundary flux 
(left ordinate) with accumulated heating (right ordinate), 
in degrees of potential temperature integrated over the FMG. 
Compensatory radiative cocling in the CMG balanced the heat- 
ing in the FMG at each vertical (pressure)level. 

meting curve merely increases monotonically, since the 
heating rate was constant in time, but the accumulated 


5۶۰۰٠٠٢٠٠۰٠٦٦٦ boundary Flux curve behaved quite differently. 


In Figure 23, portions of the flux curve with positive 
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SoBe represent flux into the FMG, while those with negative 
slope indicate outward flux. The diagonal slashes on the 
flux curve denote the two movements of the FMG. 

Two interesting features of the accumulated flux 
Seen in Figure 23 are the initial flux into the FMG, and 
EE oscillatory behavior of the flux for the remainder of 
the integration period. The initial flux into the FMG must 
be a result of adjustments between the scales in the inter- 
face region, since no ascent cell had been established at 
this early point in the integration, and the initial state 
Was barotropic. In particular, the area averages assigned 
to the internal interface coincident points, and the external 
interface values which are determined by interpolation from 
CMG tendencies, may be quite different from the initial 
values, even after a few time steps. After approximately 
eight hours, the outward flux iS monotonic, but oscillates 
with a period of about 16 hours. The reason for this oscil- 
Mti on is not clear, but it could be related to the FMG 
mewement, which occurred each 16 hours for this test. Tests 
with different diabatic heating and mean zonal flow might 
help to indicate the cause of oscillation in the thermal 
Boumdary flux. Even though the flux behaved rather erratic- 
ally, an important result was nevertheiess established. That 
the mutually determined boundary specification in the 
three-dimensional model allowed both inward and outward 


ERermal fluxes with ne Apparent computational difficulties. 
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SR TEST 2: A TEST OF INITIALIZATION AND SIMULATED 
FRONTAL DEVELOPMENT 


In the previous tropical storm case, the FMG was large 
enough to completely encompass the region of large gradients 
and heating maximum. In the second test, to be described in 
this section, however, large gradients of wind and tempera- 
ture were allowed at the interface. The objective of the 
frontal simulation test was twofold. First, a brief test of 
initialization procedures was desirable in a relatively 
smooth ,flow to perhaps suggest useful methods for real data 
initialization of the meshed system. Secondly, it was neces- 
sary for the sake of completeness of these tests of the 
meshing technique to study the behavior of the model in a 
more complex environment with a less idealized initial state. 
In this manner, some guidelines might be established as to 
size or other limitations of the 6. 

Deformation of an unstable baroclinic wave, closely re- 
sembling atmospheric frontogenesis, has been described by 
mariela (1970). This deformation resulted ےہ‎ zontd 
shear in the mean wind field when a small amplitude disturb- 
ance was Specified in the initial geopotential. This dis- 
turbance was independent of pressure, and superimposed on 
the basic zonal flow. The initial 1000 mb geopotential is 


given by 





L_ ZN 0-0 
91000 i U 37 Cos ( 7 ) sin E ) (27) 





and the potential temperature by 


R/c 
p P 
p>) A Y y Y any 28 
0 = f R ar [5 F AT SIR ( ya) (28) 
۱ 39 
F F ad al (pe ج‎ P \ 
E | i -l , -1 ad the uni- 
where o = 10 ٠۹۰ت .2م‎ OOM: 0> op > 


form lapse rate, was -0.376 K mb”, The width, or north= 


south extent, of the domain is denoted by w, and the mean 
temperature 8, was 300 K. The base pressure Py was 1000 mb 
as in previous experiments. For this test of the meshed 
model, an additional term was included in (27) which caused 
Ehe mean Zonal flow at the middle level (550 mb) to be zero. 
This was done to eliminate the necessity of movement of the 
FMG, which was of no interest in this test. Thus there was 
strong easterly flow at lower levels, replaced by westerlies 
above 550 mb. Although this pattern was somewhat unrealis- 
tic for mid-latitude flow, the dynamics were unchanged, and 
more convenient for this test of the meshed model. Other 
details of this initial state are discussed by Hadfield 


(1970), and are of no further concern here. 


Geopotential heights at the upper levels were then ob- 
tained with (23), and geostrophic winds were computed to 
complete the initialization. This simulated data was then 
the real” data for the CMG and uniform coarse grid model, 
with the 550 mb geopotential heights and temperature pattern 
depicted in Figure 24. With the FMG located as indicated by 


the dotted outline in Figure 24, two tests of initialization 
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1 6 16 21 25 
5.24. Initial 550 mb height field - 5000 m (solid) and 
potential temperature (dashed). The position of the FMG 

IT lhe initialization tests is indicated by the dotted out- 
line. 

dme performed. The first initialization test involved 
direct linear interpolation of all dependent variables from 
the CMG values shown in Figure 24 onto the FMG. No attempt 
was made to perform diagnostics of any kind, and the inte- 
dation commenced from this initial state in the FMG. Large 
amplitude, short wavelength modes were excited almost imme- 
parace ly and aliasing apparently led to non-linear computa- 
onal instability in Slichitly over 24 hours. This result 


strongly indicated that interpolated values can probably not 


be used to initialize the FMG of a multi-level meshed model, 
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the primitive equations. This is probably‏ ۲7117265 ت۰۱۲ 
particularly true when strong, nonlinear gradients of depen-‏ 
dent variables appear in the data.‏ 

The same interpolation procedure was carried out for the 
second test, except that only the winds were interpolated. 
The resulting FMG data then underwent a diagnostic phase 
after a consistent vertical motion pattern was obtained from 
(22). Balanced heights and hence temperatures were computed 
mom (24), and with no other differences from the first 
case, this integration was carried out for 24 hours. The 
resulting forecast was qualitatively very similar to that 
obtained with the uniform coarse grid. This result suggested 
that this combination of interpolation and diagnostics may 
Ké useful for real data initialization in the meshed model. 
When the initial state shown in Figure 24 is integrated for 
several model days, the small amplitude disturbance results 
"deformation of the baroclinic zone. This deformation 
Causes the formation of a frontal, or mid-latitude cyclone 
Wave, with the initially small amplitude disturbance devel- 
Oping rapidly. Very strong winds develop, and in order to 
test the computational stability aspects of the meshing 
٣۰۰٦٠٦۹ ٭‎ im a complex environment, the FMG was located as 
depicted in Figure 25. The height and potential temperature 
fields portrayed in Figure 25, are those resulting after 
five model days (600 CMG time steps) in the uniform coarse 
grid model. These predicted values were then used as data 


reo initialize the PMG. This test had two main objectives, 
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1 6 11 16 21 25 
eS 25. Results of the uniform coarse grid for a five day 
mereecast. The solid contours are 550 mb heights - 5000 m, 
and the dashed contours are isentropes. This forecast was 
the "initial" data for the meshed model test. See text for 
details. 
Ehe first of which was to subject the FMG interface to ex- 
tremely large nonlinear gradients of all dependent variables 
to determine stability characteristics and other possible 
effects. Secondly, it was desired to study the solution in 
Ee FMG after the previously described initialization tech- 
wue was utilized. In this case, however, the vertical 


Hoti on patterns were fully deyeloped, and the nonlinear terms 


ME ould include the vertical advective quantities. 
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The same initialization procedure which proved the better 
of those tested in the less perturbed model atmosphere (Fig- 
ure 24) was again employed for this test. The location of 
the FMG in Figure 25 was rather arbitrary, since not only did 
the regions of maximum potential temperature and height gra- 
dients not coincide, but the system was tilted toward the 
west at upper levels. It was obvious that the FMG depicted 
in Figure 25 was too small to model the cyclonic scale dis- 
turbance, but was nevertheless ideal oe the test described 
ES Results of a one day forecast from the "initial" 
state shown in Figure 25 are depicted in Figure 26 for the 


20 





۳۰۰ 26. Uniform Coarse grid forecast for day six at 550mb. 
the solid contours are heights - 5000 m, and the dashed 
contours are isentropes. 
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uniform coarse grid model. The 550 mb height and potential 
temperature analyses at this time show that the gradients of 
both of these variables have increased, particularly on the 
"warm front" side. Strengthening of the gradients in this 
region was also noticed by Hadfield (1970) in similar exper- 
iments with a uniform grid model. The meshed model, on the 
other hand, produced a considerably different solution, as 
portrayed in Figure 27. Although deepening of the low was 
not as large, the 5100 m contour has closed, indicating en- 
hanced development. As may be seen in Figure 27, the shape 
of both the height and potential temperature fields are dif- 
ferent for the meshed model as compared to the uniform grid 
model result in Figure 26. This is not surprising, since 
phase speed errors are much less in the FMG, and development 
ment be expected to be accelerated over that in the coarse 
id. 

Although the gradients have been strer:gthened in the FMG, 
it is difficult to argue that the solution is improved over 
e uniform grid result. At this point of the integration 
(6 days, 720 CMG time steps), the uniform coarse grid incre- 
ment is still sufficiently small to resolve the gradients, 
and possible distortion at the interface due to phase speed 
differences may mask improvements resulting from the finer 
resolution in the FMG. To be assured of a more accurate 
solution for the frontal simulation, one would need to utilize 
an FMG several times larger than that depicted in Figure 2/. 


Two important results of this test were that the mutually 
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Fig. 27. Meshed model forecast for day So mb: 

The solid lines are heights - 5000 m, and the dashed con- 
tours are isotropes. 

determined boundary specification retained the computational 


ability of the solution even in res orpable situations, 


and that phase error was reduced in the FMG. 
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vee TAE HER MODEL: A MULTIPLE GRID INCREMENT REDUCTION 


Previous sections of this report have demonstrated the 
Mutually determined boundary specification in several tests 
with a single reduction of grid scale between the coarse and 
fine grids. Although the abrupt increment change is compu- 
tationally stable, and has been shown to yield results simi- 
lar to those obtained on a uniform fine grid, one might 
argue heuristically that a gradual reduction in grid scale 
ment prove less “traumatic” to the numerical solution. 

This might be particularly true in nonlinear systems with a 
Marge grid scale reduction. In the tests of the MR model 
described below, the same basic four-level primitive equa- 


tion model discussed in previous sections was utilized. 


A. THE GRID AND MARCHING PROCESS 

Several problems are encountered with multiple grid in- 
crement reduction if the mutually deternined boundary speci- 
fication is to be employed. One important limitation is 
that the ratio of reduction must be such that certain coin- 
cident points are maintained. As may be seen in the example 
depicted in Figure 28, each successive grid increment must 
mrereasce In the Same ratio, or an integer multiple of the 
ratio, of the previous increment decrease. This is import- 
ant in the movement of the FMG, in the boundary condition 
wec cation, and dn the actual coding of the model. In 


Figure 28, the grids are labeled CMG (coarse mesh grid), 
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me.) 20. An example of a horizontal level in the MR model. 
For this example, each of the smaller grids covers 8 x 8 
increments of the next larger, and the spatial increment is 
reduced by one-half. See text for explanation of indices. 
MMG (medium mesh grid), and FMG (fine mesh grid). The indi- 
cated index points are, as in the SR model, acronyms which 
indicate the location of either the MMG or the FMG. For 
example, IWEM, indicated I, West side, ternal interface 
and Medium grid. It should perhaps be emphasized that in 
Pie grid depicted in Figure 28, the acronymic indices ending 


with M are CMG indices, and those ending with F are MMG 
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indices. Thus the presence of the FMG is not "felt" by the 
meme xcept indirectly through the MMG. Similarly, the CMG 
has no effect on the FMG other than by indirect effects 
Meemsmitted through the MMG. In this manner, the MMG acts 

as a "buffer" between the largest and the smallest scale 
modeled. Thus regions of large nonlinear gradients may be 
modeled on the small scale desired, but the expense of having 
the entire FMG very small scale is eliminated. 

The marching process for the MR model is more complex 
than the SR model, although it is analogous. Assuming that 
the first time step (which requires special attention) has 
been taken, the marching process may be described by refer- 


ring to the schematic in Figure 29. By taking one centered 


re .‏ 
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ie.) 29. A schematic to indicate the marching process in 
the MR model. The heavy line indicates any arbitrary posi- 
Zoran time after the first time step, with the time be- 
Eween each point on the line representing one FMG time 
increment. 


Pets tep, the CMC Solution may be raised from time level A 


to time level E. The tendencies around the MMG external 
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interface may then be determined by spatial and temporal in- 
En polation. A centered time step in the MMG raises the 
Solution in this grid from time level A to time level C, and 
the MMG external boundaries are updated from the previously 
m ermined CMG tendencies. Two centered time steps then 
need to be taken in the FMG to raise the solution in this 
grid from time level A to time level C. In this case, the 
external boundary values are determined from the MMG tenden- 
cies, and applied after each FMG time step. Finally, the 
internal interface values of the MMG are کت‎ dase the 
area average of the appropriate points in the FMG. The solu- 
tions in the two smaller grids are then raised from time 
Enel C to time level E by repeating the latter process. 

The final step in the marching process is to assign the CMG 
internal interface values by the appropriate area average in 
eme MMG. The procedure for the first time step is quite 
Salar, except that the first time step in each grid is 
morward in time, rather than centered. 

Movement of the refined grids in the MR model is also 
analogous to that in the SR model. As the movement initia- 
tor (e.g., surface low pressure, maximum wind gradient,etc.) 
approaches the internal interface of the FMG, the FMG moves 
eae MHC spatial increment in the appropriate direction. The 
MMG remains stationary within the CMG during this move, and 
leading and trailing edge values of the FMG at appropriate 
time levels are determined exactly as in the SR model. When 


the FMG approaches the edge of the MMG, the MMG moves one 
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MiG interval in the indicated direction, but during this move 
the FMG remains stationary. Thus the movement initiator for 
the MMG is independent of the atmospheric feature of inter- 
est, and depends only upon the location of the FMG. Leading 
and trailing edge values for the MMG are obtained in the 
usual manner after the movement, and all FMG location indices 
(i.e. IWEF, JSEF, etc.) must be redefined as the MMG "slides" 


under the stationary FMG. 


EL TROPICAL STORM SIMULATION IN THE MR MODEL 

One test similar to the tropical storm simulation in the 
SR model was performed utilizing the MR model. The major 
difference was that the mean zonal easterly flow was removed 
to eliminate the necessity of movement of the refined grids. 
The primary objectives of this test were to examine time and 
Size requirements for the different grid scales, and to 
qualitatively compare the numerical results with those ob- 
tained with the SR model, It must, of course, be emphasized 
that removal of the zonal current will affect the solution, 
but similarities between the SR and MR solutions should still 
exist. 

٦۰٦3۲۰3٦٦23۲٣٦76 for the MR model was similar to that in 
the SR case. Geostrophic winds were analytically specified 
with (25), except that F(y) = 0. The divergence equation 
(24) was then solved for the heights first in the CMG, then 
in the MMC, and finally in the FMG. Boundary values for the 


MMG and FMG are obtained in an analogous manner as in the SR 
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model. The point vortex depicted in Figure 30 was, there- 
Bere, identical in amplitude and location to the one utilized 
for the SR test, but the actual values of the heights were 
modified because the mean zonal current was removed. The 
dotted outlines in Figure 30 indicate the position of the 
refined grids, with the FMG centered in the MMG, and also 


centered over the surface low center. 





Fig. 30. Initial heights - 1400 m for the MR tropical storm 
simulation. The dotted outlines denote the location of the 
MMG (outer), and the FMG (inner). The CMG grid increment is 
indicated in the lower left corner. 

The 850 mb height forecast for 48 hours is depicted in Figure 


31 for the MR meshed system. As in the SR model, coincident 


CMG values in the MMG and FMG were replaced by those taken 
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Fig. 31. Forecast 850 mb heights - 1400 m for the MR model. 
Mie location of the refined grids is indicated by the dotted 
erreE lines. 

from the refined grids. The low center deepened approxi- 
mately 40 meters, more than double that in the SR model. 
This was probably primarily due to the removal of the basic 
mean zonal flow since all other physical parameters were 
identical. This is further demonstrated by the 700 mb ver- 
Be motion pattern in the FMG depicted in Figure 32. In 
the SR model, the vertical motion tended for form narrow 
bands of ascent, while the ascent cell in the MR migrated to 


the low center. This suggested that formation of narrow 


areas of ascent may be related to the movement of the system 
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1 5 9 13 17 


IEE. 32. Vertical motion pattern (mb x 10°) after 48 
hours at 700 mb for the FMG of the MR model. 

through the surrounding environment. The main result of the 
MR test was to demonstrate that an improved detailed solution 
could be obtained with a gradual grid scale reduction, but 
Ber etion of the computer requirements. For the meshed 
MR model depicted in Figure 28, the computer memory require- 
ment was approximately two-thirds, and the time requirement 


approximately one-half, that utilized in the p = 5 SR model 
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utilized in earlier experiments. This was in spite of the 
fact that the MR model refined grid region covered 8 x 8 CMG 


grid intervals, compared with 6 x 6 for the SR version. 
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VITIT. SUMMARY 


important EUrrent problem in meteorology is the inclu- 
sion of small scale sources and sinks of energy in numerical 
models of the atmosphere. Related to this problem is the 
question of interactions between the planetary scale flow 
and the relatively small scale features which contribute the 
energy. Design of numerical models with uniform fine grids, 
which could possibly incorporate all these scales of motion, 
is currently quite impossible, due to computer time and 
memory limitations. 'In this paper, a possible solution to 
these important problems has been proposed. Different 
spatial increments were employed in the nested system, with 
the coarser grid increment an integer multiple of the finer 
۰۶۶6ء٦۰‎ che inclusion, or nesting, of a finer resolution 
grid within a coarser one, allows small scale phenomena to 
be examined in detail without impractically large computer 
requirements. 

The initial discussion of the nesting technique dealt 
with an analytic computational stability consideration and 
the appropriate related numerical tests of a simple linear 
equation, which was solved on both uniform and nested grid 
systems. TWO unique aspects of the meshed grid system were 
the mutual interaction between the fine grid region and the 
nn COarSse grid due bo the method of specifying 


the interface boundary values, and the fact that the 
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equation was integrated simultaneously in time on both grids. 
An approach due to Platzman (1954) which utilized the ampli- 
ration matrix of the system of differential difference 
equations was used to show computational stability of the 
numerical solution when certain interface boundary specifica- 
tion was utilized. This analysis was valid when the same 
time increment was used in both grids, and the fine grid 
region remained stationary. Numerical solutions of the lin- 
ear equation obtained with the nested grid system were shown 
to compare favorably with those resulting from integration in 
a uniform fine resolution grid. The mutually determined 
boundary specification was shown to be the only satisfactory 
one of those tested when centered space and time differencing 
was used. Analytical and numerical tests of other boundary 
Seca fications showed that these specifications were either 
eemputationally unstable, or produced inferior results. 

Iwo important facets of the meshing technique also arose 
from the one-dimensional linear experiments. The first was 
Amer crucial fact that the computational stability of the 
mutually determined specification extended to the case where 
different time increments were used in each grid. This is 
particularly significant from the standpoint of obtaining a 
more accurate solution in the coarse grid, as well as most 
em icient Utilization Of computer time during the integration. 
Secondly, it was shown that movement of the fine grid region 
within the coarse grid had no apparent deleterious effect on 


the numerical solution. Subsequent numerical tests with 
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fee initial conditions, appropriate time increment in each 
grid, and movement of the fine grid to remain centered over 

a maximum gradient, again produced a numerical solution com- 
parable to that obtained on a uniform fine grid. 

Subsequent one-dimensional numerical tests with a non- 
linear system of equations demonstrated computational stabi- 
lity in more complex situations. An initial state which 
developed extremely sharp (eventually discontinuous) gradients 
was utilized to compare numerical solutions on uniform grids 
with those of the meshed system. The meshed grid solutions 
in this case were superior to those obtained on a uniform 
coarse grid, and nearly identical to those resulting when a 
form fine resolution grid was used. 

Since the ultimate goal of this research was to develop 
a meshed three-dimensional primitive equation model, exten- 
sion to a two-dimensional version was the next logical step. 
The two-dimensional model was a ten-level primitive equation 
model. The existence of internal and external gravity waves 
in such multi-ievel models was important in testing the mesh- 
Bo technique, because every interface boundary point may be 
an outflow boundary for either advective motions or gravity 
waves. The initial state for this test was a small scale 
heat source embedded in an otherwise undisturbed flow pattern, 
which stimulated the development of a convective type cell. 
Comparison of five-day forecasts showed the meshed system 
NOt Only Virtually reproduced the uniform grid solution in- 


tic Tine EIlICd, but outside this region as weil. This 
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result strongly suggested that the mutually determined speci- 
fication allowed interaction between scales at the interface 
boundaries. This was further demonstrated in favorable com- 
parisons of boundary fluxes of kinetic and thermal energy 
across the same spatial boundary in uniform fine and meshed 
systems. Certain computational devices were employed in the 
primitive equation models. The most significant was the 
lization of horizontal momentum diffusion in the fine grid 
to suppress high frequency اج تہ سس‎ noise caused by 
truncation error differences between the grids, and slight 
inconsistencies in the boundary specifications. This was 
important because such short wavelength noise may eventually 
Mead to nonlinear computational instability in primitive 
equation models. 

All three-dimensional tests were performed utilizing a 
four-level primitive equation model. The fine grid covered 
approximately ten percent of the total spatial domain in the 
three-dimensional model. Initialization of data for these 
models was achieved by means of a divergence equation, which 
was derived from the difference form of the prediction equa- 
tions. In all cases, winds were analytically specified, and 
balanced geopotential heights and temperatures were obtained 
by means of the divergence equation. The first test utilized 
an analytic initial state which simulated a tropical storm 
scale disturbance which was easily encompassed by the fine 
A region: In this Simulation, a circularly symmetrical 


vortex was embedded in the trough of a smooth wave in the 
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easterly flow. Heating was specified to roughly approximate 
the latent heat release inherent in such disturbances. Fore- 
cast results from a meshed model with a grid reduction of one- 
fifth in the fine grid were compared to those from a uniform 
coarse grid. No three-dimensional tests with uniform fine 
grids were possible because of computer limitations. The 
meshed model forecast depicted more detail of the storm simu- 
lation, and the movement of the vortex more closely agreed 
with the movement which would have resulted from purely ad- 
vective processes. This test suggested that with a carefully 
Specified initial state, and a fine grid which completely 
encompassed the region of interest, the meshed system yielded 
an improved solution. 

To test the meshed model and the initialization scheme 
in a more complex situation, the second test employed the 
same grid, but included a frontal scale disturbance too 
large to be encompassed by the fine mesh region. The initial 
data for the meshed model was a numerical solution from a 
uniform coarse grid after several days of integration, and 
thus had a fully developed vertical motion pattern. The 
data were interpolated to the fine grid, and underwent a diag- 
nostic phase with the previously derived divergence equation. 
Experiments showed that starting the integration without the 
diagnostic phase led to computational instability. An impor- 
Bantspoimtsconcernine this Initialization was that the initial 
is uerlr€e NOE NON divergent, and all nonlinear terms, includ- 


ing vertical flux, were included in the divergence equation. 
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mimes test yielded two important results. First, the mutual- 
ly determined boundary specification retained the computa- 
tional stability of the numerical solution even in cases 
where strong, nonlinear gradients occurred at the interface. 
Secondly, the solution in the fine grid of the meshed model 
after one model day beyond the initial time appeared to en- 
hance development over the entire spatial domain, and parti- 
early in the fine grid system. 

The final three-dimensional tests were performed with 
the same four-level primitive equation model, but with more 
than one reduction in grid scale. This was important be- 
cause it allowed a small region of interest to be closely 
examined with a very fine grid, but without the computer 
erme and memory expense which is inherent when the entire 
fine grid is extremely small scale. Thus, fairly active 
regions surrounding the disturbance, but some distance away, 
could still be modeled with finer resolution than the coarse 
grid, while the region of primary concern was modeled with 
even greater resolution. One test was performed with this 
meshed system, in which the initial state was a stationary 
vortex, identical to that in the single reduction tropical 
Storm simulation, except the mean easterly flow was removed. 
Qualitative comparisons of two day forecasts between the 
single reduction model and the multiple reduction model de- 
monstrated that improved solutions could be obtained with 


the multiple reduction model. These solutions compared 
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favorably with those obtained in the single reduction model, 
but utilized much less computer time and memory. 

Many possible applications of the meshing technique dis- 
cussed in this report exist. One of the most important has 
been suggested by several of the numerical experiments, 
namely forecasting the movement and development of tropical 
storms. It appears quite evident that if a fine mesh region 
was kept centered over a tropical storm in a prediction 
model, forecasts of these features would be improved. 

Another interesting related application would be experiments 
with tropical storm modification, in which the behavior of 
the storm system could be examined when heat is added at 
upper levels, for example. There are many other applications 
of meshed models, not only in meteorology, but also in other 
disciplines. When systems of equations are solved on finite 
grids, there is usually some limited portion of the solution 
domain that would be desirable to model with finer resolution 


than surrounding areas. 
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APPENDIX A 
SEMISDISCRETE AMPLIFICATION FOR MESHED 
AND UNIFORM GRID SYSTEMS 

This appendix lists the semi-discrete amplification ma- 
trices discussed in section 2a. In each case, the CMG cyc- 
lic continuity is indicated by the first and last rows of 
E matrix. The mutually determined specification is retained 
for the meshed model at the inflow boundary (fourth row of 
each meshed model matrix), and the various outflow specifi- 
cations listed in Table l are indicated at the outflow 
Boundary (tenth row). 
Matrix 1: Uniform grid with arbitrary external specifica- 


Cion. This specification is stable for all 
cases (see Platzman 1954). 
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Meshed grid, mutually determined (Specification 1 


in Table 1) 
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outflow value equal to adjacent 


uniform grid, 


Matrix 5: 


(Specification 5 in Table 1) 


interior value. 
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outílow value equal to second 


adjacent interior point value (Specification 6 


Uniform grid, 
in Table 1). 


Matrix 7: 
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Meshed grid, 


Matrix 8: 


adjacent interior point value. 
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oUt lowsevdııe determined by linear 
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extrapolation from two adjacent interior points 
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points. 
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